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Abstract 



We develop the analytic and numerical tools for data analysis of the gravitational-wave signals 
from spinning neutron stars for ground-based laser interferometric detectors. We study in detail 
the statistical properties of the optimum functional that need to be calculated in order to detect 
the gravitational-wave signal from a spinning neutron star and estimate its parameters. We derive 
formulae for false alarm and detection probabilities both for the optimal and the suboptimal filters. 
We assess the computational requirements needed to do the signal search. We compare a number of 
' criteria to build sufficiently accurate templates for our data analysis scheme. We verify the validity 

, of our concepts and formulae by means of the Monte Carlo simulations. We present algorithms by 

""^ ■ which one can estimate the parameters of the continuous signals accurately. 

PACS number(s): 95.55.Ym,04.80.Nn,95.75.Pq,97.60.Gb 

bX). 

I Introduction 

■ ■ — 1 ■ 

This paper is a continuation of the study of data analysis for one of the primary sources of gravita- 
tional waves for long-arm ground-based laser interferometers currently under construction 
spinning neutron stars. In the first paper of this series || (hereafter Paper I) we have introduced a 
two-component model of the gravitational-wave signal from a spinning neutron star and we have derived 
the data processing scheme, based on the principle of maximum likelihood, to detect the signal and esti- 
mate its parameters. In the second paper (6| (hereafter Paper II) we have studied in detail accuracies of 
estimation of the parameters achievable with the proposed data analysis method. 

The main purpose of this paper which is Paper III of the series is to study the statistical properties 
of the optimal functional that we need to calculate in order to detect the signal. We find that the two- 
component model of the signal introduced in Paper I can be generalized in a straightforward way to 
the A^-component signal. The main idea of this work is to approximate each frequency component of 
the signal by a linear signal by which we mean a signal with a constant amplitude and a phase linear 
in the parameters of the signal. We have demonstrated the validity of such an approximation in Paper 

II by means of the Monte Carlo simulations which show that the rms errors calculated using the linear 
model closely approximate those of the exact model. The key observation is that for the linear model 
the detection statistics is a homogeneous random field parametrized by the parameters of the signal. For 
such a field one can calculate a chracteristic correlation hyperellipsoid which volume is independent of the 
values of the parameters. The correlation hyperellipsoid determines an elementary cell in the parameter 
space. We find that the number of cells covering the parameter space is a key concept that allows the 
calculation of the false alarm probabilities that are needed to obtain thresholds for the optimum statistics 
in order to search for significant signals. We use these ideas to calculate the number of filters needed to 
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do the search. We show that the concept of an elementary cell is also useful in the calculation of true rms 
errors of the estimators of the parameters that can be achieved with matched filtering and explain their 
deviations from rms errors calculated from the covariance matrix. In this paper we develop a general 
theory of suboptimal filters which is necessary as such filters usually occur in practice. Our concept of an 
elementary cell carries over to the case of suboptimal filtering in a straightforward manner. The analytic 
tools develop in this work lead to independent criteria for construction of accurate templates to do the 
signal search. We demonstarte that those criteria give a consistent picture of what a suitable template 
should be. In an appendix to this paper we indicate how to parametrize the templates in order that they 
realize an approximately linear model so that the analytic formulae developed here can directly be used. 

The plan of the paper is as follows. In Sec. 2 we introduce an A-component model of the gravitational- 
wave signal from a spinning neutron star. In Sec. 3 we study in detail the detection statistics for the 
A-component model. We show that the detection statistics constitutes a certain random field. We derive 
the probabilities of the false alarm and the probabilities of detection. We present two approaches to the 
calculation of the probability of false alarm: one is based on dividing the parameter space into elemetary 
cells determined by the correlation function of the detection statistics and the other is based on the 
geometry of random fields. We compare the theoretical formulae with the Monte-Carlo simulations. In 
Sec. 4 we carry out detailed calculations of the number of cells for the all-sky and directed searches. 
In Sec. 5 we estimate the number of filters needed to calculate the detection statistics and we obtain 
the computational requirements needed to perform the searches so that the data processing speed is 
comparable to data aquisition rate. We compare our calculations with the results of Brady et al. 
obtained before by a different approach. In Sec. 6 we present in detail the theory of suboptimal filters 
and consider their use in the detection of continuous signals. In Sec. 7 we propose a detailed algorithm to 
estimate accurately the parameters of the signal and we perform the Monte-Carlo simulations to determine 
its performance. In Appendix A we give analytic formulae for some coefficients in the detection statistics. 
In Appendix B we present analytic formulae for the components of the Fisher matrix for the approximate, 
linear model of the gravitational-wave signal from a spinning neutron star. In Appendix C we give a 
worked example of the application of our theory of suboptimal filtering derived in Sec. 6. In Apendix D 
we study the transformation of the paramaters of the signal to a set of parameters such that the model 
is approximately linear. 

2 The TV-component model of the gravitational- wave signal from 
a spinning neutron star 

In Paper I we have introduced a two-component model of the gravitational-wave signal from a spinning 
neutron star. The model describes the quadrupole gravitational-wave emission from a freely precessing 
axisymmetric star. Each of the components of the model is a narrowband signal where frequency band 
of one component is centered around a frequency f a which is the sum of the spin frequency and the 
precession frequency and the frequency band of the second component is centered around 2/ D . A special 
case of the above signal consisting of one component only describes the quadrupole gravitational wave 
from a triaxial ellipsoid rotating about one of its principal axes. In this case the narrowband signal 
is centered around twice the spin frequency of the star. However there are other physical mechanisms 
generating gravitational waves and this can lead to signals consisting of many components. Recently two 
new mechanisms have been studied. One is the r— mode instability of spinning neutron stars 0, [l0| 
that yield a spectrum of gravitational- wave frequencies with the dominant one of 4/3 of the star spin. 
The other is a temperature asymmetry in the interior of the neutron star that is misaligned from the 
spin axis . This can explain that most of the rapidly accreting weakly magnetic neutron stars appear 
to be rotating at approximately the same frequency due to the balance between the angular momentum 
accreted by the star and lost to gravitational radiation. Therefore in this paper we shall introduce a 
signal consisting of A narrowband components centered around A different frequencies. More precisely 
we shall assume that over the bandwidth of each component the spectral density of the detector's noise 
is nearly constant and that the bandwidths of the components do not overlap. 

Analytic formulae in this paper will be given for the A-component signal. However in numerical 
calculations and simulations we shall restrict ourselves to a one-component model. 
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We propose the following model of the iV-component signal: 



N 4 

h{t) = Y,hi(t), hi(t) = ^A H h u {t), l = l,...,N, (1) 

1=1 i=i 

where An are AN nearly constant amplitudes. The amplitudes are nearly constant because they depend 
on the frequency of the gravitational wave which is assumed to change little over the time of observation. 
The amplitudes An depend on the physical mechanism generating gravitational waves, as well as on the 
polarization angle and the initial phase of the wave [cf. Eqs. (28)-(35) of Paper I]. The above structure 
of the TV-component signal is motivated by the form of the two-component signal considered in Paper I 
[cf. Eq. (27) of Paper I]. The time dependent functions hu have the form 



h n (t) = a(t)cos$/(i), h n (t) = b(t)cos$i(t), 
M*) = a(*) sin h H (t) = b(t)ain§i(t), 



I 1 V. (2) 



where the functions a and b are given by 

aft) = — sin27(3 — cos2A)(3 — cos2<5) cos[2(a — <fi r — H r t)] 
16 

— — cos 27 sin A(3 — cos 25) sin[2(a — 4> r — fi r f )] 
+— sin 27 sin 2A sin 25 cos [a — <\> r — fi r i] 

— — cos 27 cos A sin 25 sin [a — 4> r — fi r t] 
3 

+ - sin 27 cos 2 A cos 2 5, (3) 
b(t) — cos27sin Asin<5cos[2(a — <f> r — fl r t)] 

+— sin 27(3 — cos2A) sin5sin[2(a — <j> r — fl r t)] 
+ cos 27 cos A cos 5 cos [a — <fi r — fl r t] 

+ — sin 27 sin 2A cos 5 sin[a — <f) r — fl r t] . (4) 

The functions a and b are the amplitude modulation functions. They depend on the position of the 
source in the sky (right ascension a and declination 5 of the source), the position of the detector on the 
Earth (detector's latitude A), the angle 7 describing orientation of the detector's arms with respect to 
local geographical directions (see Sec. II A of Paper I for the definition of 7), and the phase <f> r determined 
by the position of the Earth in its diurnal motion at the beginning of observation. Thus the functions a 
and b arc independent of the physical mechanisms generating gravitational waves. Formulae (§) and (§ 
are derived in Sec. II A of Paper I. 

The phase $; of the Ith component is given by 

Ji (fc , f k+i 9 Jl,(>=)t k 9tt J^(*)t k 

= 2 - E * 7fc-w + T no ' rEs(<) £ fl B + ~ na ■ rE w X> it?' (5) 

k=0 V k=0 ' k=0 

where rgg is the vector joining the solar system barycenter (SSB) with the center of the Earth and 
joins the center of the Earth with the detector, n is the constant unit vector in the direction from 
the SSB to the neutron star. We assume that the Ith component is a narrowband signal around some 

(0) (fc) 

frequency // which we define as instantaneous frequency evaluated at the SSB at t = 0, // (k = 1, 2, . . .) 
is the A:th time derivative of the instantaneous frequency of the Zth component at the SSB evaluated at 
t = 0. To obtain formula (JsJ) we model the frequency of each component in the rest frame of the neutron 
star by a Taylor series. For the detailed derivation of the phase model see Sec. II B and Appendix A of 
Paper I. 



3 Optimal filtering for the TV-component signal 
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3.1 Maximum liklihood detection 



Maximum likelihood detection and parameter estimation method applied in Paper I to the two-component 
signal generalizes in a straightforward manner to the iV-component signal. 

We assume that the noise n in the detector is an additive, stationary, Gaussian, and zero-mean 
continuous random process. Then the data x (if the signal h is present) can be written as 

x(t) = n(t) + h(t). (6) 

The log likelihood function has the form 

\nA = (x\h)-±(h\h), (7) 
where the scalar product ( • | • ) is defined by 

{hi\h 2 ) := 4MJ — -^jj- — df. (8) 

In Eq. (|8|) ~ denotes the Fourier transform, * is complex conjugation, and Sh is the one-sided spectral 
density of the detector's noise. As by our assumption the bandwidths of the components of the signal 
are disjoint we have (hi\hi>) ~ for I ^ I', and the log likelihood ratio (Q) can be written as the sum of 
the log likelihood ratios for each individual component: 



lnA«£ 



1 

(x\h l )--{h l \h l ) 



i=i 



(9) 



Thus we can consider the iV-component signal as N independent signals. Since we assume that over the 
bandwidth of each component of the signal the spectral density Sh(f) is nearly constant and equal to 
Sh(fi), where fi is the frequency of the signal hi measured at the SSB at t — 0, the scalar products in 
Eq. (||) can be approximated by 

MM « IT^TT / ' x(tMt)dt, (MM « "rT~TTT / " Mt)} 2 dt, (10) 
&h{jl) J-T /2 bh(Jl) J-T /2 

where T is the observation time, and the observation interval is [— T /2, T a /2\. 
It is useful to introduce the following notation 
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T„/ 2 

x(t)dt, (11) 



-To/2 

Using the above notation and Eq. (|l0|) the log likelihood ratio from Eq. (^|) can be written as 

JV 

Shift) V^'" l/ 2 



^-E^rfw-I^)). d2) 



i=i 



Proceeding along the line of argument of Paper I [cf. Sec. Ill A of Paper I] we find the explicit analytic 
formulae for the maximum likelihood estimators An of the amplitudes An: 



where we have defined 



An w 
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(xhn) 
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-C(xhn) 
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[xhu) 


-C(xh m ) 
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B := 


(b 2 ), 


C := (ab) 



Z = 1,...,A, (13) 



A:=(a^), B:=(b z ), C := (ab) , D := AB — C . (14) 
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To obtain Eqs. (|13|) we have used the following approximate relations: 
(hnhi 3 ) w {hnhu) ~ (h l2 h l3 } w (h^hu) « 0, 

i = (15) 

<^l) * (A? 8 > * |A « (^4) « ^, » (^3^4) « |c, 

One can show that when the observation time T D is an integer multiple of one sidereal day the function 
C vanishes. To simplify the formulae from now on we assume that T a is an integer multiple of one sidereal 
day (in Appendix A we have given the explicit analytic expressions for the functions A and B in this case) . 
In the real data analysis for long stretches of data of the order of months such a choice of observation 
time is reasonable. Then Eqs. (|13| ) take the form 

-> (xhn) (xhn) -r (xhn) t- (xhu) 

A n ~2^-^, A l2 ^2^-^, A i3 «2^-^, A l4 w 2^-^, l = l,...,N. (16) 

The reduced log likelihood function T is the log likelihood function where amplitude parameters An 
were replaced by their estimators An. By virtue of Eqs. fll5| ) and ( |l6| ) from Eq. ( |l2| ) one gets 



2V 

2T 



(a;/i;i) 2 + (xfos) (xft, ;2 ) 2 + (xhu)" 



A B 



(17) 



To obtain the maximum likelihood estimators of the parameters of the signal one first finds the 
maximum of the functional T with respect to the frequency, the spindown parameters and the angles a 
and 6 and then one calculates the estimators of the amplitudes An from the analytic formulae (|l3| ) with 
the correlations (xhu) evaluated at the values of the parameters obtained by the maximization of the 
functional T . Thus filtering for the TV-component narrowband gravitational-wave signal from a neutron 
star requires AN linear filters. The amplitudes An of the signal depend on the physical mechanisms 
generating gravitational waves. If we know these mechanisms and consequently we know the dependence 
of An on a number of parameters we can estimate these parameters from the estimators of the amplitudes 
by least-squares method. We shall consider this problem in a future paper. 

Next we shall study the statistical properties of the functional T . The probability density functions 
(pdfs) of T when the signal is absent or present can be obtained in a similar manner as in Sec. Ill B of 
Paper I for the two-component signal. 

(fc> 

Let us suppose that filters hu are known functions of time, i.e. the phase parameters fi, a, 8 are 
known, and let us define the following random variables: 

x u :={xhu), l = l,...,N, i = l,...,4. (18) 

Since x is a Gaussian random process the random variables xu being linear in x are also Gaussian. Let 
Eo{a;;i} and Eijxji} be respectively the means of xu when the signal is absent and when the signal is 
present. One easily gets 

E o {%} = 0, i = l,...,4 ) l = l,...,N, (19) 

Ex{z u } = §^L4 U , Ei{xu} = %BA 12 , E 1 {x l3 } = ±AA l3 , Ei{x u } = \BA U , l = l,...,N. (20) 

Since here we assume that the observation time is an integer multiple of one sidereal day it immediately 
follows from Eqs. ( |T5| ) that the Gaussian random variables xu are uncorrelated and their variances are 
given by 

V«r{*„}=^^ ) i=l,3, 

stnn 1 = 1,...,*. (21) 

Var M = ^p, * = 2,4, 

The variances are the same irrespectively whether the signal is absent or present. We introduce new 
rescaled variables z;, ; : 



zu = 2 




1 = 1,. ..,N, (22) 
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so that zu have a unit variance. By means of Eqs. (|l^) and ( pp[ ) it is easy to show that 

E {z u } = 0, i = l,.. .,4, Z = l,...,/V, (23) 



and 




l = l,...,N. (24) 



ToB_ 
Shift) 



-An, 

The statistics .F from Eq. ( |l7| ) can be expressed in terms of the variables zu as 

iV 4 



=i i=i 



The pdfs of T both when the signal is absent and present are known. When the signal is absent 2!F 
has a x 2 distribution with 4N degrees of freedom and when the signal is present it has a noncentral % 2 
distribution with 47V degrees of freedom and noncentrality parameter A = Y]r—-i 2t=i TO H' ^ e nn< ^ that 
the noncentrality parameter is exactly equal to the optimal signal-to-noise ratio d defined as 



d := y/(h\h). (26) 
This is the maximum signal-to-noise ratio that can be achieved for a signal in additive noise with the 



linear filter |12 . This fact does not depend on the statistics of the noise. 



Consequently the pdfs po and p\ when respectively the signal is absent and present are given by 

Po(F) = [y exp(-n (27) 
r 2 jr\(n/2-l)/2 . / _ 1 



Pi(d,F) = ] dnl2 _ x I n/2 - 1 (dV2T)eK V (-T--d 2 y (28) 

where n — AN is the number of degrees of freedom of \ 2 distributions and I n /2~i is the modified Bessel 
function of the first kind and order n/2 — 1. The false alarm probability Pp is the probability that T 
exceeds a certain threshold T Q when there is no signal. In our case we have 

P F {T )~ / PG {T)dT = eM-Fo) IT' ( 29 ) 

k=0 

The probability of detection Pp, is the probability that T exceeds the threshold T when the signal-to- 
noise ratio is equal to d: 

/>oo 

P D {d,T ) := J Pl {d,T)dT. (30) 

The integral in the above formula cannot be evaluated in terms of known special functions. We see 
that when the noise in the detector is Gaussian and the phase parameters are known the probability of 
detection of the signal depends on a single quantity: the optimal signal-to-noise ratio d. 

Our signal detection problem is posed as the statistical hypothesis testing problem. The null hypothesis 
is that the signal is absent from the data and the alternative hypothesis is that the signal is present. The 
test statistics is the functional T . We choose a certain significance level a which in the theory of signal 
detection is the false alarm probability defined above. We then calculate the test statistics T and compare 
it with the threshold T calculated from equation a = Pp(JF ). If J- exceeds the threshold T we say 



G 



that we reject the null hypothesis at the significance level a. The quantity 1 — a is called the confidence 
level. Clearly because of the statistical nature of the problem the null hypothesis can be rejected even if 
the signal is present. In the theory of hypothesis testing we call the false alarm probability the error of 
type I and the 1 — Pjj which is the probability of false dismissal of the signal we call the error of type II. 
When the signal is known by Neyman-Pearson lemma the likelihood ratio test is the most powerful test 
i.e. it maximizes the probability of detection Pjj which in the theory of hypothesis testing is called the 
power of the test. 



3.2 False alarm probability 

Our next step is to study the statistical properties of the functional T when the parameters of the phase of 
the signal are unknown. We shall first consider the case when the signal is absent in the data stream. Let 
£ be the vector consisting of all phase parameters. Then the statistics T(£) given by Eq. (|l7]) is a certain 
generalized multiparameter random process called the random field. If the vector £ is one-dimensional the 
random field is simply a random process. A comprehensive study of the properties of the random fields 
can be found in the monograph [ p"3| . For random fields we can define the mean m and the autocovariance 
function C just in the same way as we define such functions for random processes: 

m(£) := E{JF(£)}, (31) 

C(£,£') := E{TO)-m(€)][^(0-m(0]}- (32) 

We say that the random field T is homogeneous if its mean m is constant and the autocovariance function 
C depends only on the difference £ — £ . The homogeneous random fields defined above are also called 
second order or wide-sense homogeneous fields. 

In a statistical signal search we need to calculate the false alarm probability i.e. the probability that 
our statistics T crosses a given threshold if the signal is absent in the data. In Paper I for the case of a 
homogeneous field T we proposed the following approach. We divide the space of the phase parameters £ 
into elementary cells which size is determined by the volume of the characteristic correlation hypersurface 
of the random field T . The correlation hypersurface is defined by the requirement that the correlation 
C equals half of the maximum value of C . Assuming that C attains its maximum value when £ — £ = 
the equation of the the characteristic correlation hypersurface reads 

C(r) = ic(0), (33) 

where we have introduced r :=£ — £'. Let us expand the left hand side of Eq. ( |33| ) around r = up to 
terms of second order in r. We arrive at the equation 

M 

J2 G tJ r lTj = 1, (34) 



where M is the dimension of the parameter space and the matrix G is defined as follows 

1 d 2 C{r) 



ij : ~ C(0) ()-<)-, 



(35) 



The above equation defines an M -dimensional hyperellipsoid which we take as an approximation to the 
characteristic correlation hypersurface of our random field and we call the correlation hyperellipsoid. The 
M-dimensional Euclidean volume V^ e ll of the hyperellipsoid defined by Eq. (f34j) equals 

n M/2 

Kell = ; , (36) 

r(M/2 + l)\/detG V ' 

where Y denotes the Gamma function. 

We estimate the number N c of elementary cells by dividing the total Euclidean volume Vtotal of the 
parameter space by the volume V^ e ii of the correlation hyperellipsoid, i.e. we have 

Nc = J^taJ_ (37) 

Vcell 
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We approximate the probability distribution of J-(£) in each cell by probability po{!F) when the 
parameters are known [in our case by probability given by Eq. The values of the statistics T in 

each cell can be considered as independent random variables. The probability that T does not exceed 
the threshold T in a given cell is 1 — Pf{T^), where Pp(To) is given by Eq. (p9|). Consequently the 
probability that T does not exceed the threshold T Q in all the N c cells is [1 — Pp(!F )} Na . The probability 
Pp that T exceeds T Q in one or more cell is thus given by 

pT{T ) = l-[l-P F {To)] Na . (38) 

This is the false alarm probability when the phase parameters are unknown. The expected number of 
false alarms Np is given by 

Nf = N c P f {To)- (39) 
By means of Eqs. (|2^) and ([}?]), Eq. ([59|) can be written as 

cc " k=0 



Using Eq. fl39|) we can express the false alarm probability Pf from Eq. ( p8|) in terms of the expected 
number of false alarms. Using lim n _ >00 (l + — )™ = exp(x) we have that for large number of cells 

Pj(JP )«l-exp(-iV F ). (41) 

When the expected number of false alarms is small (much less than 1) we have Pf ~ Np. 

Another approach to calculate the false alarm probability can be found in the monograph jlij . Namely 
one can use the theory of level crossing by random processes. A classic exposition of this theory for the 
case of a random process, i.e. for a one-dimensional random field, can be found in Ref. fl5| . The case 
of M-dimensional random fields is treated in ||l3f and important recent contributions are contained in 
Ref. |p^| . For a random process n(t) it is clear how to define an upcrossing of the level u. We say 
that n has an upcrossing of u at t D if there exists e > such that n(t) < u in the interval (t a — e,t Q ), 
and n{t) > u in (t ,t + e). Then under suitable regularity conditions of the random process involving 
differentiability of the process and the existence of its appropriate moments one can calculate the mean 
number of upcrossings per unit parameter interval (in the one-dimensional case the parameter is usally 
the time t and n(t) is a time series). 

For the case of an A/-dimensional random field the situation is more complicated. We need to count 
somehow the number of times a random field crosses a fixed hypersurface. Let !F(£) be Af-dimcnsional 
homogeneous real- valued random field where parameters £ = (£i, ■ ■ ■ , £m) belong to Af -dimensional 
Euclidean space M. M and let C be a compact subset of R M . We define the excursion set of inside 
C above the level T as 

Ar{T , C):={CeC: > T } . (42) 

It was found that when the excursion set does not intersect the boundary of the set C then a suitable 
analogue of the mean number of level crossings is the expectation value of the Euler characteristic \ °f 
the set Ajr. For simplicity we shall denote xl-Ari^o, C)] by xt„- It turns out that using the Morse 
theory the expectation value of the Euler characteristic of A? can be given in terms of certain multidi- 



mensional integrals (see Ref. 13 , Theorem 5.2.1). Closed form formulae were obtained for homogeneous 
M-dimensional Gaussian fields and 2-dimensional x 2 fields (see |T^], Theorems 5.3.1 and 7.1.2). Recently 
Worsley jl6| obtained explicit formulae for M-dimensional homogeneous x 2 field. We quote here the 
most general results and give a few special cases. 

We say that [/(£), £ e M M , is a X 2 field if U{£) = £f=i A ; (|) 2 , where X x (£),... , X n (£) are inde- 
pendent, identically distributed, homogeneous, real-valued Gaussian random fields with zero mean and 
unit variance. We say that E/(£) is a generalized x 2 field if the Gaussian fields Xi(£) are not necessarily 
independent. 

Let 2J r (^) be a x 2 field and let Xi(£), I = 1, . . . , n, be the component Gaussian fields then under 
suitable regularity conditions (differentiability of the random fields and the existence of appropriate 
moments of their distributions) 



nxr ] = M^ft, ^"' M)/2 exp{-Fo)W M>n (r o ). (43) 
tt-m/^I (n/2) 



In Eq. (43) V is the volume of the set C and matrix A is denned by 



a _ <9 2 C(0 



(44) 





= 1, 








= T 


-\{n- 


-1), 


w 3>n 


= n 


-(n- 


\)T + ±(n-l)(n-2), 




= n 


+ |(n 


- 1) 2 ^ 2 - |n^ - i(n - l)(n - 2)(n - 3) 



where C is the correlation function of each Gaussian field Xi(£). Wu,n(T ) is a polynomial of degree 
M — 1 in T given by 

v ; i=o fc=o v 

where division by factorial of a negative integer is treated as multiplication by zero and [N] denotes the 
greatest integer < N. We have the following special cases: 



(46) 



It has rigorously been shown that for the homogeneous Gaussian random fields the probability distri- 
bution of the Eulcr characteristic of the excursion set asymptotically approaches a Poisson distribution 
(see Ref. |Q, Theorem 6.9.3). It has been argued that the same holds for x 2 fields. It has also been shown 
for M-dimensional homogeneous x 2 fields that asymptotically the level surfaces of the local maxima of 
the field are M-dimesional ellipsoids. Thus for large threshold the excursion set consists of disjoint and 
simply connected (i.e. without holes) sets. Remembering that we assume that the excursion set does not 
intersect the boundary of the parameter set the Euler characteristic of the excursion set is simply the 
number of connected components of the excursion set. Thus we can expect that for a x 2 random field 
the expected number of level crossings by the field i.e. in the language of signal detection theory the 
expected number of false alarms has a Poisson distribution. Thus the probability that f max does not 
cross a threshold T is given by exp(— E[x^J) and the probability that there is at least one level crossing 
(i.e. for our signal detection problem the false alarm probability Pp) is given by 

P^{T Q ) = PCFma* > To) w 1 - cxp(-E[x^J). (47) 

From Eqs. (^l|) and ([|7]) we see that to compare the two approaches presented above it is enough to 
compare the expected number of false alarms Nf with E[x^ D ]. It is not difficult to see that for x 2 fields 
G = 2A. Thus asymptotically (i.e. for large thresholds To) using Eqs. (|6|), (f4(]), and ( ff3"l ) we get 



V . 2 M / 2 r(M/2 + l)T- M ' 2 as T Q -» oo, (48) 



E[x^, 



where we have used that V from Eq. (|43|) coincides with Vtotai from Eq. (|40|) . 

Worsley ([^6), Corollary 3.6) also gives asymptotic (i.e. for threshold T a tending to infinity) formula 
for the probability P(T mayL > T ) that the global maximum T max of T crosses a threshold T Q '- 



P(T max > To) -> J t i^ 2) ^ n+M)/ ^^M-To) as To^oc. (49) 

In the signal detection theory the above probability is simply the false alarm probability and it should 
be compared with the probability given by Eq. (^). It is not difficult to verify that asymptotically the 
Eqs. ( p8| ) and (^) are equivalent if we replace expected number of false alarms Np by E[xjfJ. This 
reinforces the argument leading to Eq. d4§|). 

The above formulae were obtained for continuous stationary random fields. In practice we shall always 
deal with a discrete time series of finite duration. Therefore to see how useful the above formulae are in 
the real data analysis of discrete time series it is appropriate to perform the Monte Carlo simulations. 
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We have first tested Eqs. (|27]) and (|2^) for the probability density of the false alarm and the false 
alarm probability in the simplest case of n = 2 and the known signal. Using a computer pseudo-random 
generator we have obtained a signal x consisting of TV = 2 s independent random values drawn from a 
Gaussian distribution with zero mean and unit variance. The optimal statistics J- in this case is 

Fk = ^j-, k=l,...,N/2+l, (50) 

where \Xk\ is the modulus of the fcth component of the discrete Fourier transform of x. In other words the 
optimal statistics is the periodogram sampled at Fourier bins. When x consists of independent identically 
distributed Gaussian random variables we know jl7j that for k = 2, . . . , N/2 the statistics 2Tk has a \ 2 
distribution with 2 degrees of freedom whereas for k = 1 (zero frequency bin) and k — N/2 + \ (maximum, 
Nyquist frequency bin) has a x 2 distribution with 1 degree of freedom. In our Monte Carlo simulation 
we have generated the signal 10 6 times and we have made histograms of 127 bins of the statistics Th- In 
the upper plot of Figure 1 we have shown the (appropriately normalized) histogram for all the Fourier 
bins for k — 2, . . . , N/2 and in the middle plot of Figure 1 we have presented the cumulative distribution. 
Thus the probability density generated in the upper plot is to be compared with Eq. ( p7| ) for n = 2 
whereas the distribution obtained in the middle one is to be compared with Eq. (^9|) for n = 2. Both 
theoretical distributions are exponential and they are given by solid curves in the two plots. The last bin 
in the upper plot of Figure 1 deviates substantially from the exponential curve. This is because in this 
bin all the events above the maximum value of the histogram range are accumulated. We get 11 events 
altogether in the last bin. The expected value of the events calculated from Eq. ( |39| ) is 8.25 (where we 
have put N c — 127). In the two lower plots of Figure 1 we have presented the cumulative distributions for 
the first and the last bin. The theoretical cumulative distribution that follows from the \ 2 distribution 
with 1 degree of freedom is given by 1 — erf( 1 /^ r /2) (solid curve in the plots). We see that simulated 
and theoretical distributions agree very well. 

We have next tested the formulae for the false alarm probabilities given by Eqs. ( |38| ) and (^) against 
the Monte Carlo simulations. We have considered again the case of n — 2 and we have simulated the 
optimal statistics for the case of a monochromatic signal (M — 1) and the case of a signal with one 
spindown included (M — 2). We have generated the random sequence of length N = 2 s as in the 
first simulation described above. We have however introduced an extra parameter P — the zero padding. 
Namely we add zeros to the random sequence so that its total length is (1 + P)N. When we take Fourier 
transform of the zero-padded signal we get additional points in the Fourier domain between the Fourier 
bins. Zero padding essentially amounts to interpolating the periodogram between the Fourier bins. Thus 
the larger the P the closer the discretelly sampled periodogram to a continuous function. To generate 
the statistics T for the signal with one spindown we have multiplied the generated random sequence 
x(k) by T(k; I) = exp[-27riZ<7i(fc - l) 2 ], where k = 1, . . . , N ,1 = 10, . . . , 29, and a x = (3/tt) ^/oT^. The 
function T(k; I) is called a filter or a template. The multiplication operation is the matched filtering 
which in our case is also called dechirping. The quantity o\ is the accuracy of estimation of the 1st 
spindown parameter for the optimal signal-to-noise ratio d = 1 divided by \[2 and it is the maximum 
extent of the ellipse defined by Eq. (|34|) measured from the origin along the spindown axis in the Cartesian 
(frequency, spindown)-planc. The parameter o\ defines the spacing of the templates that we choose in 
our simulations. The zero padding is always done after the dechirping operation. Our optimal statistics 
is the modulus of the discrete Fourier transform of the dechirped and zero padded data divided by the 
number of points in the original data (2 8 in our case). We have made 10 5 trials. 

The results of the simulation are presented in Figure 2. The three upper plots are the results for the 
monochromatic signal search and the three lower ones are the results for the one spindown signal search. 
The false alarm curves are given in the plots on the left. We see that the false alarm probability exhibits 
a threshold phenomenon, it drops very sharply within a narrow range of the detection threshold. We see 
from the plots on the left of Figure 2 that for P = (no zero padding) the results of the simulation agree 
well with Eq. ( |3S| ) (solid line) whereas for P = 3 there is a reasonable agreement with Eq. (|4^) (dashed 
line) . 

In the plots on the right of Figure 2 we have divided the probability of the false alarm obtained 
from the simulations by the probabilities obtained from the theoretical formulae. The upper plots give 
comparison with Eq. (p8| ) based on dividing the parameter space into cells whereas the lower plots give 
comparison with Eq. (E7I) based on the expectation value of the Euler characteristic of the excursion set. 
We see that for the monochromatic signal for thresholds up to 10 Eq. (|3|) gives a reasonable agreement 
for P = whereas Eq. ( fl?]) gives a good agreement for P = 3. For the frequency modulated signal for 



10 



P = Eq. ([38|) underestimates the false alarm probability whereas Eq. ((47|) overestimates the false alarm 
probability. For P = 3 there is an underestimate of the false alarm probability by both formulae. For 
thresholds greater than 10 the curves become irregular what may be attributed to a sparse number of 
events for such large thresholds. 



3.3 Detection probability 

When the signal is present a precise calculation of the pdf of T is very difficult because the presence of the 
signal makes the data random process x(t) nonstationary. As a first approximation we can estimate the 
probability of detection of the signal when the parameters are unknown by the probability of detection 
when the parameters of the signal are known [given by Eq. (|30|)1. This approximation assumes that when 
the signal is present the true values of the phase parameters fall within the cell where T has a maximum. 
This approximation will be the better the higher the signal-to-noise ratio d. Parametric plot of probability 
of detection vs. probability of false alarm with optimal signal-to-noise ratio d as a parameter is called the 
receiver operating characteristic (ROC). 

We have performed the numerical simulations to see how the ROC obtained from the analytical 
formulae presented above compares with that obtained from the discrete finite duration time series. We 
have generated the noise as in the simulations of the false alarm probability and we have added the 
signal. We have considered both the monochromatic and the linearly frequency modulated signal. The 
frequency of the signal was chosen not to coincide with one of the Fourier frequencies. However in the 
dechirping operation to detect the frequency modulated signal we have chosen the spindown parameter 
in the filter to coincide with the spindown parameter of the signal. We have perfomed 10 trials and we 
have examined the cumulative distributions of the two Fourier bins between which the true value of the 
frequency of the signal had been chosen. The results are presented in Figure 3. The two upper plots are 
for the monochromatic signal and the lower two are for the one spindown signal. In the plots on the left 
we compare the probability of detection calculated from Eq. (^) with the results of the simulations and 
in the plots on the right we compare the theoretical and the simulated receiver operating characteristics. 
For the false alarm probability we have used the formula (|3Sj) . In the inserts we have zoomed the ROC 
for small values of the false alarm probability. We see that the agreement between the theoretical and 
simulated ROC is quite good. 



4 Number of cells for the one- component signal 

Let us return to the case of a gravitational-wave signal from a spinning neutron star. In Sec. 5 of Paper 
II we have shown that each component of the ./V-component signal can be approximated by the following 
one-component signal: 

h{t; h a , $o,l) = K sin [$(t; £) + $ ] , (51) 
where the phase $ of the signal is given by 

= y^fc ( Tfr ) "I i a i [Res sin ((f> + fl a t) + .Re cos A cos e sin (<£ r + Q r t)\ 

k=o VJo ^ c 

+a 2 [Res cos (0 O + ft a t) + Re cos A cos (</> r + Sl r t)]} , (52) 



where T a denotes the observation time, Res — 1 AU is the mean distance from the Earth's center to the 
SSB, Re is the mean radius of the Earth, fl Q is the mean orbital angular velocity of the Earth, and <f) 
is a deterministic phase which defines the position of the Earth in its orbital motion at t = 0, e is the 
angle between ecliptic and the Earth's equator. The vector £ collects all the phase parameters, it equals 
£ = (ai, <X2, u>q, . . . , <x> s ), so the phase $ depends on s + 3 parameters. The dimcnsionlcss parameters 

w _ 
are related to the spindown coefficients f a introduced in Eq. (m as follows: 

The parameters ct\ and et 2 arc defined by 

a\ := f (cos e sin a cos S + sin e sin 8) , a% := f cos a cos S. (54) 



11 



In Appendix D we show that the parameters ax, «2 can be used instead of the parameters a, 6 to label 
the templates needed to do the matched filtering. 

The signal defined by Eqs. ( f3l| ) and ([52]) has two important properties: it has a constant amplitude 
and its phase is a linear function of the parameters In Paper II we have shown that for this signal's 
model the rms errors calculated from the inverse of the Fisher information matrix reproduce well the rms 
errors of the full model presented in Sec. 2. We will use here the simpified signal (|5l])-(|5^) to estimate 
the number of elementary cells in the parameter space. 

For the signal given by Eqs. ( pl| ) and ( |52| ) the statistics T of Eq. (jlT]) can be written as 

^(0^{M0] 2 + M€)] 2 }> (55) 

where 



a?c(€):=2W-^y (a; cos $(*;€)), x s {£) := 2 J sin . (56) 



We calculate the autocovariance function C [defined by Eq. (|32|)1 of the random field ( |55| ) when the 
data x consists only of the noise n. We recall that n is a zero mean stationary Gaussian random process. 
Consequently we have the following useful formulae jl8| 

E{(n|hi)(n|h a )} = (fti|ft 2 ), (57) 

E{(n|h 1 )(n|h 2 )(n|h 3 )(n|h 4 )} = (hi\h 2 )(h 3 \h 4 ) + (hi\h 3 )(h 2 \h 4 ) + (hi\h 4 )(h 2 \h 3 ), (58) 
where hi, h 2 , h 3 , and h 4 are deterministic functions. Let us also observe that 

(sin 2 $(i;0)«(cos 2 $(;;0)«^ (59) 



Making use of Eqs. (57), (|5q), and (59) one finds that 

E {^(€)} « 1, (60) 
E « 1 + 2 [(cos $(t; £) cos $(f ; £')} 2 + (cos $(t; £) sin $(t; O}' 

+ (sin *(t; |) cos $(t; O}* + (sin *(t; £) sin £')) 2 ] , (61) 

where subscript means that there is no signal in the data. For our narrowband signal to a good 
approximation we have 

(cos$(t;0cos$(t;O) « 5 (cos[*(t; - $(t; £')]} , (62) 

(cos sin $(t;£')) « ~\ (sin[$(t;0 - *(*;0]> > (63) 

<Bin*(t;€)coB*(t;0) ~ 5 <™[*(t;€) - <*>(*;£')]> , (64) 

(sin$(i;0sin$(t;O) « ^ <oos[*(t; €) - $(t; OJ> ■ (65) 



Collecting Eqs. (J6QD — (J65[) together one gets 

C(t,f) = Eo{^(OT£')}-Eo{^(€)}Eo{^(0} 

« (cos[$(t;0-^(i;O]> 2 + <sin[$(t;C)-^(*;O]> 2 - (66) 

The phase $ given by Eq. (^32j) is a linear function of the parameters £ hence the autocovariance 
function C from Eq. ( |66[ ) depends only on the difference r = £ — £' and it can be written as 

C(r) « (cos[$(t; r)]) 2 + <sin[$(i; r)]) 2 . (67) 

To calculate the volume of the elementary cell by means of Eq. ( |36| ) we need to compute the matrix 
G defined in Eq. @. Substituting @ into @ we obtain 

G = 2f, (68) 
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where the matrix T has the components 

The matrix T is the reduced Fisher information matrix for our signal where the initial phase parameter 
$o [cf. Eq. (pi])] has been reduced, see Appendix B. 

As the mean ( |60| ) of the random field T is constant and its autocovariance ( |66| ) depends only on the 
difference £ — £' the random field is a homogeneous random field. Let us observe that for the fixed 
values of the parameters £ the random variables x c and x s are zero mean and unit variance Gaussian 
random variables. However the correlation between the Gaussian fields x c and x s does not vanish: 

E{x c (£)x s (0} * (sin[$(t;0 - S(t;£')]) , (70) 

and thus the Gaussian random fields x c and x s are not independent. Therefore T is not a \ 2 random field 
but only a generalized \ 2 random field. Our formula for the number of cells [Eq. (pi|)] and the formula 
for the false alarm probability [Eq. d3Sj)] apply to any homogeneous random fields however formula ( p3[ ) 
applies only to \ 2 fields. Nevertheless by examining the proof of formula (^) [|l3| [l6| we find that it is 
very likely that the formula holds for generalized x 2 random fields as well if we replace the determinant 
of the matrix A by the determinant of the reduced Fisher matrix T. 

The total volume of the parameter space depends on the range of the individual parameters. Following 
Ref. [Q we assume that 

27rT / min < luq < 27rT / max , (71) 
1 / T Q \ k 

-/3fe^o < u k < Pk^o, where (3 k ■= , , : — — , k = 1, . . . , s, (72) 



k + 1 \t,____ , 

where / m i n and / max are respectively the minimum and the maximum value of the gravitational-wave 
frequency, T m ; n is the minimum spindown age of the neutron star. The parameters a± and a 2 defined in 
Eq. (|54|) fill, for the fixed value of the frequency parameter loq, 2-dimensional ball concentrated around 
the origin in the (ai, a2)-plane, with radius equal to loq / '(2nT ): 



(ai,aa)GB 2 (0,wo/(27iT )). (73) 

-all 

total(s) 



Taking Eqs. (|7l|)-(|73|) into account the total volume K 3,1 . 1 , of the parameter space for all-sky searches 



V ^L(,) = I duj I I da x da 2 I du) X ... / duj s 



with s spindowns included can be calculated as follows 

""total(s) ~~ 

J J J B 2 (0,LJo/(2irT )) 

27rT / min -fcuo -(3 s u> 

= — T s+1 — ) (f+ 3 -ff a ) (74) 

(« + 3)(«+l)! ° \T min J Umax JmmJ [ > 

The volume V^, a \ of one cell we calculate from Eq. (^) for M = s + 3 and G = 2r?R , where iVl is 
the reduced Fisher matrix ( |69| ) for the phase $ given by Eq. (^2|) with s spindowns included: 

KeV) = ^ , (75) 



r(( S + 5)/2)^detr^ 

In Appendix B we have given formulae needed to calculate matrices T^' for s = 0, . . . , 4 analytically. In 
Figure [| we have plotted the volume V^L^. of one cell as a function of the observation time for signals 
with various numbers s of spindown parameters included. 

The number A^Jw s ) of cells for all-sky searches is given by 

Aral] _ total(s) _ A 7T / „ i c _ s+1 / , s+3 fS+3 N , , 

coUs(s) ""^7~ r (s /2 + i) V detr (*) {^-) T ° (/— ( 76 ) 
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In Figure |5| we have plotted the number AT^} lg , , of cells as a function of the observation time T for various 
models of the signal depending on the minimum spindown age r m i n and the maximum gravitational- 
wave frequency / max , and for various numbers s of spindowns included (the minimum gravitational- wave 
frequency / m i n = 0). We see that for a given T min and / max curves corresponding to different numbers 
s intersect. This effect was observed and explained by Brady et al. Q. To obtain the number of cells 
for a given observation time T we always take the number of cells given by the uppermost curve. We 
have calculated the observation times T^ ss , k -, for which the numbers of cells with k and k + 1 spindowns 
included coincide: 

-^cc!ls(fc+l) {To = ^cross(fc)) = ^cc!ls(fc)(^o = T CTOSa ^), k = 0, . . . , S — 1. (77) 

In Table | we have given the values of T^,^ for all the signal models considered. 



T m i„ (years) 


/max (Hz) 


T c1L (k) ( da ys) 


TZss (k) (days) 


k = 


k = 1 


k = 2 


k = 3 


k = 


k = 1 


k = 2 


fe = 3 


40 


10 3 


0.21 


3.11 


116 


311 


0.03 


3.53 


40.5 


175 


40 


200 


0.31 


5.19 


158 


389 


0.06 


6.04 


60.5 


242 


10 3 


10 3 


0.46 


114 


575 


2210 


0.14 


30.2 


452 


2300 


10 3 


200 


0.69 


157 


725 


3040 


0.32 


51.7 


676 


3180 



Table 1: The observation times for which the numbers of cells with k and k + 1 spindowns included 
coincide for various models of the signal depending on the minimum spindown age T m ; n and the maximum 
gravitational-wave frequency / max . The minimum gravitational-wave frequency / m i n = 0. In the case of 
all-sky searches we have used the latitude A = 46.45° of the LIGO Hanford detector and we have put 
(f> = 0.123 and <p r = 1.456. 

The Fisher matrix depends on the phases 4> r , (j> , and the latitude A of the detector (see Appendix 
B). We know from Paper II (see Appendix C of Paper II) that the Fisher matrix also depends on the 
choice of the instant of time at which the instantancnous frequency and spindown parameters are defined 
(in the present paper this moment is chosen to coincide with the middle of the observational interval). We 
find that the determinant detT^ and consequently the number of cells does not depend on this choice. 
The dependence on the remaining parameters is studied in Figure ^. The dependence on the phases <j) r 
and 4> is quite weak. The dependence on A is quite strong however for the detectors under construction 
for which A varies from 35.68° (TAMA300) to 52.25° (GEO600) the number of cells changes by a factor 
of 2 for 7 days of observation time and by around 10% for 120 days of observation time. 

In Sec. 5 of Paper II we have shown that for directed searches the constant amplitude signal given by 
Eqs. ( |5"l| ) and ([52]) can be further simplified by discarding in the phase ([32]) terms due to the motion of 
the detector w.r.t. the SSB and the rms errors calculated form the inverse of the Fisher matrix do not 
change substantially. Such a signal reads 



s / t \ k+1 

h(t;h o ,$ ,£) = ft„ Bin [*(*;£) + *„], *(*;£) = l> fc (jf- ) ( 78 ) 



The vector £ has now s + 1 components: £ = (luq, . . . , lu s ). 

Using Eqs. (|Fl]) and ( [72] ) the total volume V^jj^x of the parameter space with s spindowns included 
for directed searches is calculated as follows 

27rTo/ max 010)0 P.wo 

Kotal( S ) = J du) Q J dwi... J du s 

2irTofmin -0iu) O -Ps^o 

s(s+l)/2 



(»+!)(* + 1)1 



\ 'min / 



The volume V t d ^ al{s) of one cell we calculate from Eq. @ for M = s + 1 and G = 2T^, where rgj 
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is the reduced Fisher matrix (| 



for the polynomial phase (|78|) with s spindowns included: 
.. (V2) (s+1)/2 



V"' 



dl(s 



r(( s + 3)/2)^detr^ 



(80) 



The matrix for s = 0, . . . , 4 can be calculated analytically by means of formulae given in Appendix 



B. 



The number A^^w^) °f independent cells is given by 



ydir 

Tydir _ total(s) 
cells(s) — T/dir 

"ceU(a) 



2(35+1)72^8/2+1 

(s + l)r(s/2+ 1) 



detrg 



T 



s(s+l)/2 



J o \J max J min / 



(81) 



In Figure |^ we have plotted the number of cells N^, s ^ as a function of the observation time T a for various 
models of the signal depending on the minimum spindown age T mm , the maximum gravitational- wave 
frequency /max, and the number s of spindowns included (the minimum gravitational- wave frequency 
/min = 0). We see that like for all-sky searches for a given T nl - m and / max curves corresponding to 
different numbers s intersect. We have calculated analytically the observation times T^ se ^ for which 
the numbers of cells with k and k + 1 spindowns included coincide: 



Ardir 
Jv cclls(M 



rpdir 
cross (fc) 



/vrdir 
iV cclls(Ai) 



(T. 



cross(/c) 



0,...,s-l. 



(82) 



Using Eq. (|81|) one obtains 



cross(fc) 



(* + 2)r((* + 3)/2) 



2V2^(fc + l)r((fc + 2)/2) ^ detr^ +1) / m +x 2 - /, 



det Tfl\ 
(fe) 



ffc+i 

' max 



f 



fe+1 



mm fc+1 
^•fe+2 min 
' min 



l/(fc+2) 



1. (83) 



In Table | we have given the values of T^" s , k , for all the signal models considered. 

In Table 2 we have given the number of cells both for all-sky and directed searches for various models of 
the signal depending on the minimum spindown age T mm and the maximum gravitational- wave frequency 
/max, and for the observation time T Q of 7 and 120 days (the minimum gravitational- wave frequency 
/min = 0). The number of cells is calculated from Eq. (|7^) for all-sky searches and from Eq. ( |8l| ) in 
the case of directed searches. For a given observation time T Q the number s of spindowns one should 
include in the signal's model is obtained as such number s chosen out of s = 0, . . . , 4 for which AT?^,, 
(or A^[ ls ( s is the greatest. 

We have also calculated the threshold T for the 1% false alarm probability (or equivalently for 99% 
detection confidence). By means of Eqs. ( p9|) and (|3^) for n = 2 (what corresponds to a one-component 
signal) the relation between the threshold T Q and the false alarm probability a reads 



T = -\n \l - (1 - a) l > N ° 



a = 0.01, 



(84) 



where N c is the number of cells. Following the relation between the expectation value of the optimum 
statistics when the signal is present and the signal-to-noise ratio which is given by 



E 1 {T} = l + U 2 , 



(85) 



we have calculated the "threshold" signal-to-noise ratio 



d Q := ^2{T - I), 



(86) 



where T is given by Eq. (Q) . The values of d D for various models of the signal and observation times 
of 7 and 120 days are given in Table 2. If the signal-to-noise ratio is d Q then there is roughly a 50% 
probability that the optimum statistic will cross the threshold T a - 
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'mm 


J max 


all-sky 


directed 


s 


cells ( s) 


d 


s 


ceUs(sj 


d 


7 


40 


10 3 


2 


4.2 x 10 18 


9.6 


2 


3.7 x 10 14 


8.6 


7 


40 


200 


2 


1.3 x 10 15 


8.8 


2 


2.9 x 10" 


8.0 


7 


10 3 


10 3 


1 


1.5 x 10 ly 


9.0 


1 


1.9 x 10" 


8.0 


7 


10 3 


200 


1 


2.4 x 10 i3 


8.3 


1 


7.6 x 10 iU 


7.6 


120 


40 


10 3 


3 


3.8 x 10 2y 


12 


3 


7.2 x 10 23 


11 


120 


40 


200 


2 


1.1 x 10 2B 


11 


3 


1.2 x 10 21 


10 


120 


10 3 


10 3 


2 


2.2 x 10 2b 


11 


2 


6.0 x 10 iY 


9.4 


120 


10 3 


200 


1 


2.7 x 10 22 


11 


2 


4.8 x 10 ib 


8.9 



Tabic 2: Number of cells for all-sky and directed searches for various models of the signal depending 
on the minimum spindown age r m i n and the maximum gravitational- wave frequency / max , and for the 
observation time T of 7 and 120 days. The minimum gravitational-wave frequency / m i n = 0. To calculate 
the Fisher matrix Tfy we have used the latitude A = 46.45° of the LIGO Hanford detector and we have 
put 4> = 0.123 and 4> r = 1.456. For each case we also give the 99% confidence threshold signal-to- noise 
ratio d calculated by virtue of Eq. (^) . 

5 Number of filters for the one-component signal 

To calculate the number of FFTs to do the search we first need to calculate the volume of the elementary 
cell in the subspace of the parameter space defined by ujq = const. This subspace of the parameter space 
is called the filter space. 

Let us expand the autocovariance function C of Eq. fl67| ) around r = up to terms of second order 
in t : 

M 

C(r) = l- ^Tur^j, (87) 



where are defined in Eq. (|69|) and M is the number of phase parameters. In Eq. (|87[) we have used 
the property that C attains its maximum value of 1 for t = 0. Let us assume that t\ corresponds to 
frequency parameter and let us maximize C given by Eq. (^) with respect to t\ . It is easy to show that 
C attains its maximum value, keeping T2, ■ ■ ■ , tm fixed, for 



M 

_ 1 
Tl 



1 M 



Let us define 



C(t 2 ,...,t m ) :=C(t 1 ,t 2 ,...,t m ). (89) 
Substituting Eqs. @ and (§|) into Eq. (||) we obtain 

M 

C{t 2 , . . . , t m ) = 1 - Yj FijnTj, (90) 

where _ _ 

Ly := Tij ~ -. (91) 
I n 

We define an elementary cell in the filter space by the requirement that at the boundary of the cell 
the correlation C equals 1/2: 

C(T 2 ,...,r M ) = i. (92) 

Substituting ( po|) into ( |92| ) we arrive at the equation describing the surface of the elementary hyperellipsoid 
in the filter space: 

M 



i,j=2 
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The volume of the elementary cell is thus equal to [cf. Eq. ( p6| 

F cell = ^ _. (94) 

r((M+l)/2)VdetT 

The volume V cc ii of the elementary cell in the filter space is independent on the value of the frequency 
parameter. 

Taking Eqs. (|7l|)-(|73|) into account the total volume I^totai(s) °f the filter space for all-sky searches 
with s spindowns included can be calculated as follows 




B 2 (0,Wo/(25rT )) J J 



^totai( s ) = ill daida.2 / dui... / du s 



(95) 

w =27rT o / max 
. . a ll 

Putting in Eq. (^5|) u>o — 27rT / max we have defined V total ( s ) as that slice uiq — const of the parameter 

space which has maximum volume. 
a ii 

The volume V ce ii( s ) of one cell in the filter space we calculate from Eq. (94) for M = s + 3: 



-all _ (7r / 2 )(s+2)/2 

r(( s + 4)/2) v /detrS 



VZ Ks) = , ,„ , 07) 



where the matrix I\ s ) is calulated from Eq. (H|) forf = rgL 
The number A/|jJ eI , s , s N of filters for all-sky searches is given by 



aU _ ^totalW _ 2 (S + 2) / -all / J Q \ s+2 

JV filters( s ) - ^li - / + jypTT^ + 7^ V det L M \ Z ) MM 
* ccll(s) 



In Figure [§] we have plotted the number ^V^iters(s) °^ filters as a function of the observation time T Q for var- 
ious models of the signal depending on the minimum spindown age r m i n and the maximum gravitational- 
wave frequency / max , and for various numbers s of spindowns included. We see that for a given r m i n 
and /max curves corresponding to different numbers s intersect. This effect was observed and explained 
in Ref . : in the regime where adding an extra parameter reduces the number of filters the parameter 
space in the extra dimension extends less than the width of the elementary cell in this dimension. To 

obtain the number of filters for a given observation time T D we always take the number of filters given by 

a ii 

the uppermost curve. We have also calculated the observation times T CIOSS ^ for which the numbers of 
filters with k and k + 1 spindowns included coincide: 

^fi!tors(fc+l) (^o = 7 1 C ross(fc)) = ^filters(fc) = ^cross(fc) ) 7 = 0, . . . , S — 1. (99) 

_ a ll 

In Table || we have given the values of T cross ^ for all the signal models considered. 

— dir 

For directed searches the total volume V iotll u s \ of the filter space with s spindowns included we 
calculate using Eqs. fln]) and @: 




I3 b uj q 

^totaiW = I J dux... J duj s 



2 2s - K s ( T ^ s ( s + 1 )/ 2 



(s + 1)! 

w =27rT o / max 



(/maxT ) S . (100) 



The volume VJmi^ of one cell in the filter space for directed searches with s spindowns included we 



calculate from Eq. 1941) for M = s + 1 



F- (s) = ^ =, (ioi) 

r(( s + 2)/2)^detr ( ;; 
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T m in (years) 


/max \ nz ) 


^cross(fe) (days) 


^cross(fc) (days) 


fc = 


k = 1 


k = 2 


k = 3 


fc = 1 


fc = 2 


fc = 3 


40 


10 3 


0.19 


3.01 


113 


307 


3.26 


38.8 


171 


40 


200 


0.30 


5.07 


156 


384 


5.58 


58.0 


236 


10 3 


10 a 


0.45 


111 


566 


2170 


27.9 


434 


2250 


10 3 


200 


0.66 


153 


715 


2990 


47.7 


649 


3100 



Table 3: The observation times for which the numbers of filters with k and k + 1 spindowns included 
coincide for various models of the signal depending on the minimum spindown age T m ; n and the maximum 
gravitational-wave frequency / max . In the case of all-sky searches we have used the latitude A = 46.45° 
of the LIGO Hanford detector and we have put <j) Q = 0.123 and <fr r = 1.456. 



where the matrix T^ is calulated from Eq. 



lh for r 



■pdir 

1 w 



The number A/|jJ erB , s x of filters in the case of directed searches is thus given by: 



Ai-dir 

iv filters(s) 



V 



total(s) 



V. 



■dir 
cell(s) 



2 (3 S -2)/2 7r (s+l)/2 

r((a + 3)/2) 



detr 



s(s+l)/2 



(/» 



(102) 



In Figure ^| we have plotted the number of filters for various models of the signal depending on the min- 
imum spindown age r m j n and the maximum gravitational- wave frequency / max , and for various numbers 

— dir 

s of spindowns included. We have also calculated analytically the observation times T cross ( fc ) for which 
the numbers of filters with fc and k + 1 spindowns included coincide: 



iv filtcrs(fe+l) 



T n = T r . 



s(k) 



= N, 



filtcrs(fc) 



To = T c 



s(fc) 



fc = 1, . . . , S — 1. 



(103) 



Using Eq. (102) one obtains 



cross(/c) 



r((fc + 4)/2) 



2 v / 2lfr((fc + 3)/2)\ detrf +1) /. 



det T (fe) 



_k+l 



l/(fc+2) 



(104) 



In Table || we have given the values of T cross ^ for all the signal models considered. 

In Table 4 we have given the number of filters both for all-sky and directed searches for various 
models of the signal depending on the minimum spindown age T m - ln and the maximum gravitational-wave 
frequency / ma x, and for the observation time T D of 7 and 120 days. The number of filters is calculated from 
Eq. ( |98| ) for all-sky searches and from Eq. (102) in the case of directed searches. For a given observation 
time T a the number s of spindowns one should include in the signal's model is obtained as such number 
s chosen out of s = 0, . . . , 4 for which NS 



( or N mcrs(s)) is the greatest. 



rail 

v filters(s) 

We shall next compare the number of filters obtained above with the number of filters calculated by 
Brady et al. . In their calculations they have assumed a constant amplitude of the signal however they 
have used a full model of the phase. To calculate the number of templates they have used so called metric 
approach of Owen They have assumed a certain geometry of spacing of the templates: combination 
of a hexagonal and a hypercubic spacing and they have introduced an additional parameter — a mismatch 

which was the measure of the correlation of the two neighbouring templates. Also in their calculation 
they have assumed that the data processing method involves resampling of the time series so that the 
rcsamplcd signal is monochromatic. We shall compare the number of filters in Table 4 of our paper with 
the corresponding number of filters given in Table 1 of |Q . Our calculations correspond to mismatch 
/i = 0.5. This means that to compare our numbers of filters with the corresponding numbers of Brady et 
al. our numbers have to be multiplied by 2.4, 5.8, 15, and 40 for the signal with 0, 1, 2, and 3 spindowns 
respectively for all-sky searches and by 1.3, 1.7, 2.2 for 1, 2, and 3 spindowns respectively for directed 
searches. The difference in the volume of our hyperellipsoidal cells and their volumes of elementary 
patches means [see Ref. Q], Eq. (5.18) for all-sky searches and the paragraph above Eq. (7.2) for directed 
searches] that our numbers aditionally have to be multiplied by 1.7, 2.2, 2.8, and 3.6 for all-sky searches 
and by 1.0, 1.4, 1.3 for directed searches for comparison. After introducing the corrections for the 
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T ■ ( VPS) T"Q ) 


J max \ L1Zj J 


all-sky 


directed 


s 


Tyall 

filtcrs(s) 


V (Tf) 


s 


/ycLir 
filtcrs(s) 


V (Tf) 


7 


40 


10 3 


2 


1.4 x 10 1U 


2.6 x 10 3 


2 


9.5 x 10 5 


1.7 x 10" 1 


7 


40 


200 


2 


2.3 x 10 7 


7.8 x 10- 1 


2 


3.8 x 10 4 


1.3 x 10~ 3 


7 


10 3 


10 3 


1 


4.6 x 10 7 


8.5 


1 


3.8 x 10 3 


7.0 x 10~ 4 


7 


10 3 


200 


1 


3.7 x 10 b 


1.3 x 10~ 2 


1 


7.7 x 10^ 


2.6 x 10~ b 


120 


40 


10 3 


3 


8.4 x 10 iy 


1.7 x 10 13 


3 


1.3 x 10 i4 


2.7 x 10 Y 


120 


40 


200 


2 


1.1 x 10 1Y 


4.3 x 10 y 


3 


1.0 x lO" 


3.9 x 10 4 


120 


10 3 


10 3 


2 


4.4 x 10 lb 


9.2 x 10 8 


2 


9.0 x 10 Y 


1.8 x 10 


120 


10 3 


200 


1 


2.4 x 10 13 


9.3 x 10 b 


2 


3.6 x 10 b 


1.4 x 10- 1 



Table 4: Number of niters for all-sky and directed searches for various models of the signal depending 
on the minimum spindown age T m i n and the maximum gravitational- wave frequency / max , and for the 



—all 



observation time T of 7 and 120 days. To calculate the Fisher matrix ]?(**) we have used the latitude 
A = 46.45° of the LIGO Hanford detector and we have put <j) = 0.123 and <f> r — 1.456. For each case 
we also give the number V of flo ating point operations per second (flops) needed to do the search; V is 
calculated by means of Eq. (105). 



mismatch parameter and the size of an elementary cell we find that our corrected number of templates 
is greater than the number of templates given in Table 1 of by (going from top to bottom of Table 1) 
2.8 x 10 4 , 14, 2.7, and 1.5 for all-sky searches and by 2.2, 1.7, 0.31, and 0.25 for directed searches. We 
thus conclude that considering the differences in the way the calculations were done there is a reasonable 
agreement between the number of filters obtained by the two approaches except for one case: all-sky 
searches with the maximum frequency / max = 200 Hz and the minimum spindown age r m j n = 1000 years 
where the difference is 4 orders of magnitude. 

We would also like to point out to the uncertainties in the calculation of the number of filters. Our 
model of the intrinsic spin frequency evolution of the neutron star is extremely simple: we approximate 
the frequency evolution by a Taylor series. In reality the frequency evolution will be determined by 
complex physical processes. The size of the parameter space is likewise uncertain. The range for the 
spindown parameters [see Eqs. (f7l"|)-(f72"|)] was chosen so that the total size of our parameter space is the 
same as in M. The approximation of the time derivative of the frequency as /max/imin that is used to 
estimate the maximum value of the spindowns is probably an order of magnitude estimate. This implies 
that the size of the parameter space and consequently the number of filters is accurate within s(s + l)/2 
orders of magnitude, where s is the number of spindowns in the phase of the signal. Even this large 
uncertainty does not change the conclusion that all-sky searches for 120 days of observation time are 
computationally too prohibitive. 

To estimate the computational requirement to do the signal search we adopt a simple formula [see 
Eq. (6.11) of JtJ] for the number V of floating point operations per second (flops) required assuming that 
the data processing rate should be comparable to the data aquisition rate (it is assumed that fast Fourier 
transform (FFT) algorithm is used): 

V = 6/ max iV / [log 2 (2/ max r o ) + 1/2], (105) 

where Nf is the number of filters. The above formula assumes that we calculate only one modulus of the 
Fourier transform. Calculation of the optimal statistics T for the amplitude modulated signal requires 
two such muduli for each component of the signal [see Eq. (99) of Paper I, we assume that the observation 
time is an integer multiple of the sidereal day so that C — 0] and several multiplications. Moreover if 
dechirping operations are used instead of resampling, the data processing would involve complex FFTs. 
All these operation will not increase the complexity of the analysis i.e. the number of floating point 
operations will still go as 0(N\og 2 (N)), where N is the number of points to be processed. 

In Table 4 we have given the computer power V (in Terafiops, Tf) required for all the cases considered. 
We see that for 120 days of observation time all-sky searches are computationally too prohibitive whereas 
for directed searches only one case (r m j n = 1000 years, / m a X = 200 Hz) is within reach of a 1 Terafiops 
computer. For 7 days of observation time all cases except for the most demanding all-sky search with 
Tmin = 40 years and / max = 1 kHz are within a reach of a 1 Terafiops computer. 
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Finally we would like to point to a technique that can distribute the data processing into several 
smaller computers. We shall call this technique signal splitting. We can divide the available bandwidth 
of the detector (/ m i n , /max) into M adjacent intervals of length B. We then apply a standard technique 
of heterodyning. For each of the chosen bands we multiply our data time series by exp(— 2nifi), where 
fx = /min + IB (J = 0, . . . , M — 1). Such an operation moves the spectrum of the data towards zero by 
frequency //. We then apply a low pass filter with a cutoff frequency B and we resample the resulting 
sequence with the frequency 2B. The result is M time series sampled at frequency 2B instead of one 
sampled at 2/ max . The resampled sequencies are shorter than the original ones by a factor of M and 
each can be processed by a separate computer. We only need to perform the signal splitting operation 
once before the signal search. The splitting operation can also be performed continuously when the data 
are collected so that there is no delay in the analysis. The signal splitting does not lead to a substancial 
saving in the total computational power but yields shorter data sequencies for the analysis. For example 
for the case of 7 days of observation time and sampling rate of 1 kHz the data itself would occupy around 
10 GB of memory (assuming double precision) which is available on expensive supercomputers whereas 
if we split the data into a bandwidth of 50 Hz so that sampling frequency is only 100 Hz each sequence 
will occupy 0.5 GB memory which is available on inexpensive personal computers. 

In the case of a narrowband detector, e.g. for the GEO600 detector tuned to a certain frequency /„ 
around a bandwidth B, it is natural to apply the above data reduction technique so that the resulting 
sampling frequency is 2B. Such a technique is applied in data preprocessing of bar detectors [p3| . From 
the formulae given in the present section one can show that to perform the all-sky search with the 
integration time of 7 days for pulsars in the bandwidth of 50 Hz around the frequency of 300 Hz and 
minimum spindown age of 40 years so that the processing proceeds at the rate of data aquisition requires 
a 1 Teraflops computer. Since the data sequence occupies only 0.5 GB of memory the data processing 
task can be distributed over several smaller computers. If we also relax the requirement of data processing 
to be done in real time the signal search can be performed by a 20 Gigaflops workstation in a year. 



6 Suboptimal filtering 

It will very often be the case that the filter we use to extract the signal from the noise is not optimal. 
This may be the case when we do not know the exact form of the signal (this is almost always the case 
in practice) or we choose a suboptimal filter to reduce the computational cost and simplify the analysis. 
We shall consider here an important special case of a suboptimal filter that may be usful in the analysis 
of gravitational- wave signals from a spinning neutron star. 



6.1 General theory 

We shall assume a constant amplitude one-component model of the signal. Then the optimal (maximum 
likelihood) statistics is given by Eq. ( |55| ) . Let us suppose that we do not model the phase accurately and 
instead of the two optimal filters cos [3>(t; £)] and sin [3>(t; £)] we use filters with a phase 3>'(i; £'), where 
function $' is different form $ and the set of filter parameters £' is in general different from £, i.e. J- su b 
has the form [cf. Eqs. (|55|) and (Eg|)] 



sub 



2T 

S h (fo) 



(xcos$'(t;0) 2 + (a:sin$'(i;0) 2 > ( 106 ) 



where we have assumed that the suboptimal filters are narrowband at some "carrier" frequency f a as in 
the case of optimal filters. 

Let us first establish the probability density functions of J- su b when the phase parameters £ are 
known. Since the dependence on the data random process is the same as in the optimal case the false 
alarm and detection probability densities will be the same as for the optimal case i.e. 2J- SU ^ has a central 
or a noncentral x 2 distribution with 2 degrees of freedom depending on whether the signal is absent or 
present. From the narrowband property of the suboptimal filter we get the following expressions for the 
expectation values and the variances of J- S ub (0 means that signal is absent and 1 means that signal is 
present): 

E {^sub} = 1, Ej {T suh } = 1 + ^d 2 uh , (107) 
Var {JF sub } = 1, Van {jF sub } = 1 + d 2 uh , (108) 
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where 

1 /2 

d suh := d{(cos[<&(t;£) - $'(i;r)]) 2 + <sin[*(t;€) - $'(^')]) 2 } , (109) 

here d is the optimal signal-to-noise ratio. 

We see that for the suboptimal filter introduced above the false alarm probability has exactly the same 
X 2 distribution as in the optimal case whereas the probability of detection has noncentral \ 2 distribution 
but with a different noncentrality parameter d 2 uh . We shall call d su b (the square root of the noncentrality 
parameter) the suboptimal signal-to-noise ratio. It is clear that when the phases of the signal and the 
suboptimal filter are different the suboptimal signal-to-noise ratio is strictly less and the probability of 
detection is less than for the optimal filter. 

When the parameters are unknown the functional J-^ub is a random field and we can obtain the 
false alarm probabilities as in the case of an optimal filter. Here we only quote the formula based on 
the number of independent cells of the random field. One thing we must remember is that the number 
of cells for the suboptimal and the optimal filters will in general be different because they may have a 
different functional dependence and a different number of parameters. Thus we have [cf. Eqs. ( p9| ) and 
(||) for n = 2] 

Pj F {Fo) = 1 - [1 - exp(-T )} N ^, (110) 

where N sc is the number of cells for the suboptimal filter. 

The detection probability for the suboptimal filter is given by [cf. Eqs. (p8h and (||(]) for n = 2] 



PsD{d su h,Fo) ■= I p s l(d S ub, -T 7 ) cL? 7 , (111) 

where 



p s i{d snh ,F) = J (d snh V2T^j cxp - i 



-F-^dU)- (H2) 

Probability of detection for the suboptimal filter is obtained from the probability of detection for the 
optimal one by replacing the optimal signal-to-noise ratio d by the suboptimal one d su b- 

When we design a suboptimal filtering scheme we would like to know what is the expected number 
of false alarms with such a scheme and what is the expected number of detections. As in the optimal 
case the expected number N s p of false alarms with suboptimal filter is given by [cf. Eqs. ( p9| ) and ( |39"1 ) 
for n = 2] 

N sF = N sc eiq>{-F ). (113) 

To obtain the expected number of detections we assume that the signal-to-noise ratio d varies inversely 
proportionally to the distance from the source and that the sources are uniformly distributed in space. 
We also assume that the space is Euclidean. Let us denote by d\ the signal-to-noise ratio for which 
the number of events is one. Then the number of events corresponding to the signal-to-noise ratio d is 
(di/d) 3 . The expected number of the detected events is given by 



N D {di,F )=Z j x 2 P D (jjr^o) dx 



(114) 



in the case of the optimal filter, and by 



N sD (d lsuh ,T ) =3^ ^p^^l^L,^ dx (115) 



for the suboptimal filter. Let us note that [cf. Eq. (109)] 



1/2 

disub - di {(cos[$(f;€) - I^O]) 2 + (sin[$(i;0 - $'(<;0]) 2 } • (116) 

Because of the statistical nature of the detection any signal can only be detected with a certain 
probability less than 1. In the case of Gaussian noise for signals with the signal-to-noise ratio around 
the threshold this probability is roughly 1/2 and it increases exponentially with increasing signal-to-noise 
ratio. In Appendix C we give a worked example of the application of the statistical formulae for the 
suboptimal filtering derived above. 
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6.2 Fitting factor 

To study the quality of suboptimal filters (or search templates as they are sometimes called) one of the 
present authors |2^, |2l| introduced an I— factor defined as the square root of the correlation between the 
signal and the suboptimal filter. It turned out that a more general and more natural quantity is the 
fitting factor introduced by Apostolatos 22 . The fitting factor FF between a signal h(t; 9) and a filter 
h'(t; 6') (6 and 9' are the parameters of the signal and the filter, respectively) is defined as 

bb := max — — — (H'j 

91 ^{h{t-e)\h{t-e))^{h>{t-o')\h>(t-e')) 

If both the sig nal h and the filter h! are narrowband around the same frequency f a the scalar products 
(•]•) from Eq. (117) can be computed from the formula 



[h s h-j) * I h x (t)h 2 (t)dt, (118) 

-T /2 



S h (fo) 



where Sh is the one-sided noise spectral density and T a is the observation time. 
Let us assume that the signal and the filter can be written as 

h(t;d) = h sm^(t;C), h'{t\ 6') = h' sin$'(t; £'), ( 119 ) 

where h and h' a are constant amplitudes, £ and C' denote the para meters ent ering the phases $ and 
of the signal and the filter, respectively. We substitute Eqs. (11E) into Eq. ( 117 ). Using Eq. ( 118 ) we 
obtain ^ 

FF re max -L [ cos [$(t; C) - C')] dt - (120) 

It is easy to maximize the FF ( 120 ) with respect to the initial phase of the filter. Let us denote the initial 
phases of the functions ^ and by <&o and & , respectively. Then 

¥(t;C) = *(*;€) + *o, * , (t;0 = *'(*;0 + *o, ( 121 ) 

where £ and £ denote the remaining parameters of the signal and the filter, respectivley. After substitu- 
tion Eqs. (121) into Eq. (120) we easily get 

FF « max(cos[$(t;C)-*'(*;0 + (*o-*o)]) 

= max {cos ($ - *{,) (cos [$(t; £) - *'(t; £')] ) - sin ($ - *o) < sin C) - €')]>} 

= max {(cos [d>(t; - $'(t; £')] >' + (sin [<&(*; - *'(t; £')] • (122) 

Thus we obtain that the FF is nothing else but the ratio of the maxi mize d value of the suboptimal 
signal-to-noise ratio d su b and the optimal signal-to- noise ratio d [cf. Eq. (|l09| )1 . We stress however that 
the value of the fitting factor by itself is not adequate for determining the quality of a particular search 
template — one also needs the underlying probability distributions (both the false alarm and the detection) 
derived in the previous subsection. This is clearly shown by an example in Appendix C. 

In the remaining part of this subsection we shall propose a way of approximate computation of the 
fitting factor. Let us now assume that the filter and the signal coincide, i.e. <f>' = $, and the filter 
parameters £' differ from the parameters £ of the signal by small quantities A£: £ = £ + A£. The Eq. 
(122) can be rewritten as 

1 /2 

FF ps max { (cos [3(t; £) - $(t; £ + A£)}) 2 + (sin [$(*; £) - $(t; £ + A£)]) 2 } . (123) 



Obviously the FF (123) attains its maximu m val ue of 1 when A£ = 0. Let us expand the expression in 
curly brackets on the right-hand side of Eq. ( |123|) w.r.t. A£ around A£ = up to terms of second order 
in A£. The result is 

FF « J 1 - mm ( £ r y ■ J | , (124) 
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where 



d<f> d<f> 



a$ \ / 9$ 



(125) 



One can employ the formula (124) to estimate the FF in the case when the filter is obtained from 
the signal $ by replacing some of the signal parameters by zeros, provided the signal <E> depends weakly 
on these discarded parameters. Let the signal <E> depend on n parameters .. ,£„, and the filter $' is 
defined by 

*'(*;£,...,&) :=<K*; £!,•••, 4^^, (126) 

n — k 

where k < n, so the filter <£' depends on k parameters . . . One can write 

- *'(*; O - *(*; 6, • • • , £») - & • • • , &, o, • . . , q) = *(*; - $(*; £ + A£) (127) 



i—k 



with 



A& = 



£-6 



i — 1, . . . , fc, 
j = fc + 1, . . . , n. 



(128) 



We want to approximate the differnce $(i; £) — $(t; £ + A£) with A£ given by Eq. (128) by its Taylor 
expansion around A£ = 0. It is reasonable provided the two following conditions are satisfied. Firstly, the 
filter parameters differ slightly from the respective parameters of the signal, i.e. the quantitites A£j are 
small compared to for i = 1, . . . , k. Secondly, the function <f> depends on the parameters £fc+i, . . . , £ n 
(discarded from the filter) weakly enough to make a reasonable approximation by Taylor expansion up 
to A£,; = — £j for i = k + 1, . . . , n. If the above holds, one can use the formula (124) to approximate the 



FF. Taking Eqs. fll27|) and (128) into account, from Eq. (124) one gets 



FF 



mm 

A«i,...,A£ fc 



A{fc + i — — £fc-t 




1/2 



(129) 



6.3 Fitting factor vs. 1/4 of a cycle criterion 

Let us consider the phase of the gravitational-wave signal of the form [cf. Eq. (hh] 



4Kw t k+1 
$(i) = 2vr^/ 



S3 , <. k H k 



k=0 



(fc + 1)! 



— n • r ES (i) ^ /ory + — n ■ r E (t) ^ /„-. 



fe=0 



fc=0 



(130) 



In Paper I we have introduced the following criterion: we exclude an effect from the model of the signal in 
the case when it contributes less than 1/4 of a cycle to the phase of the signal during the observation time. 
In Paper II we have shown that if we restrict to observation times T a < 120 days, frequencies f Q < 1000 



Hz, and spindown ages r > 40 years, the phase model (13C) meets the criterion for an appropriate choice 
of the numbers si, S2, and S3. We have also shown that the effect of the star proper motion in the phase 
is negligible if we assume that the star moves w.r.t. the SSB not faster than 10 3 km/s and its distance to 
the Earth r a > 1 kpc. In Table ^|, which is Table 1 of Paper II, one can find the numbers s\, S2, and S3 
needed to meet 1/4 of a cycle criterion for different observation times T , maximum values / ma x of the 
gravitational-wave frequency, and minimum values r m i n of the neutron star spindown age. 

In Appendix A of Paper I we have indicated that the 1/4 of a cycle criterion is only a sufficient 
condition to exclude a parameter from the phase of the signal but not necessary. In this subsection 
we study the effect of ne glect ing certain parameters in the template by calculating FFs. We employ 
the approximate formula (12£) developed in the previous subsection t o ca lculate FF between the one- 
component constant amplitude signals with the phases given by Eq. ( |l30|) for numbers s\, 82, and S3 
taken from Table [5] and the same signals with a smaller number (as compared to that given in Table |s|) 
of spindowns included. We have found that for the first two models of Table [| if in the template one 
neglects the fourth spindown, FF is greater than 0.99, both for all-sky and directed searches. For other 
cases m Table I we have found that neglecting any spindown parameter can result in the FF appreciably 
less than one. 
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T (days) 


Tmin (years) 


/max (Hz) 


8 1 


S2 


S3 


120 


40 


10 3 


4 


3 





120 


40 


200 


4 


2 





120 


10 3 


10 3 


2 


1 





120 


10 3 


200 


2 


1 





7 


40 


10 3 


2 


1 





7 


40 


200 


2 


1 





7 


10 3 


10 3 


1 


1 





7 


10 3 


200 


1 


1 






Table 5: The number of spindown terms needed in various contribution s to the phase of the signal 
depending on the type of population of neutron star s searched for [cf. Eq. (|l30| )1. The number s\ refers 
to the dominant polynomial in time term in Eq. ( 13C ) , S2 refers to the Earth orbital motion contribution, 
and S3 refers to the Earth diurnal motion contribution. 



In Paper II we have considered the effect of the proper motion of the neutron star on the phase of 
the signal assuming that it moves uniformly with respect to the SSB reference frame. We have found 
that for the observation time T Q = 120 days and the extreme case of a neutron star at a distance r a 
■10 pc moving with the transverse velocity |v ns ^| = 10 3 km/s (where v ns j_ is the component of the star's 
velocity v ns perpendicular to the vector no), gravitational- wave frequency f a — 1 kHz, and spindown age 
r = 40 years, proper motion contributes only ~4 cycles to the phase of the signal. We have shown in 
Paper II that in this extreme case the phase model consistent with the 1/4 of a cycle criterion reads [cf . 
Eq. (33) in Paper II] 

- 27r E^°7FTr7! + - n ° ' rEs(<) E f°u + ~ ( n ° ' rE(<) + ' rEs(<) 1 ) U (131) 



(k + l)\ c " '^'1! c 

fe=0 v ' fe=0 



The ratio v ns j_ /r Q determines the proper motion of the star and can be expressed in terms of the proper 
motions fi a and p,$ in right ascension a and declination S, respectively (see Sec. 4 of Paper II). 



For the extreme case described above we have applied formula (129) t o ca lculate the FF between the 
one-component constant amplitude signal with the phase given by Eq. (131) and the same signal with 
a simplified phase. We have found that when both proper motion parameteres fi a , [i$ and the fourth 

(4) 

spindown parameter f a are neglected, the FF is greater than 0.99 for both all-sky and directed searches. 
Thus we conclude that neglecting the fourth spindown and the proper motion does not reduce appreciably 
the probability of detection of the signal. 

It is also interesting to compare the results obtained from the calculation of the fitting factor with 
the results summarized in Table 1 for the observation times when the number of cells for models with k 
and k + l spindowns coincides. The observation times given in Table 1 can be interpreted as observation 
times at which one should include the k+l parameter in the template. We see that for the first two 
models in Table || the Table 1 says that only 3 spindowns are needed as indicated by the calculation of 
the FF. The remaining cases also agree except for the cases of 120 days of observation time and 200 Hz 
frequency where Table 1 indicates one less spindown than Table |^. Finally we note that the crossover 
observation times in Table 1 agree within a few percent with those for the number of filters given in Table 
3. 



7 Monte Carlo simulations and the Cramer-Rao bound 

As signal-to-noise ratio goes to infinity the maximum-liklihood estimators become unbiased and their 
rms errors tend to the errors calculated from the covariance matrix. The rms errors calculated from the 
covariance matrix are the smallest error achievable for unbiased estimators and they give what is called 
the Cramer-Rao bound. 

In this section we shall study some practical aspects of detecting phase modulated and multiparameter 
signals in noise and estimating their parameters. For simplicity we consider the polynomial phase signal 
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with a constant amplitude. Our aim is to estimate the parameters of the signal accurately. We compare 
the results of the Monte Carlo simulations with the Cramer-Rao bound. 

We consider a monochromatic signal and signals with 1, 2, and 3 spindown parameters. In our 
simulations we add white noise to the signals and we repeat our simulations for several values of the 
optimal signal-to-noise ratio d. To detect the signal and estimate its parameters we calculate the optimal 
statistics T derived in Sec. 3. The maximum likelihood detection involves finding the global maximum 
of T . Our algorithm consists of two parts: a coarse search and a fine search. The coarse search involves 
calculation of T on an appropriate grid in parameter space and finding the maximum value of T on the 
grid and the values of the parameters of the signal that give the maximum. This gives coarse estimators 
of the parameters. Fine search involves finding the maximum of T using optimization routines with the 
starting value determined from the coarse estimates of the parameters. The grid for the coarse search 
is determined by the region of convergence of the optimization routine used in the fine search. We have 
determined the regions of convergence of our optimization routines in the noise free case. For the case of 
a monochromatic signal when T depends only on one parameter (frequency) our optimization algorithm 
is based on golden section search and parabolic interpolation. For a signal with some spindowns included 
T depends on s + 1 parameters (s is number of spindowns) and we use Nelder-Mead simplex algorithm. 

To perform our simulations we have used MATLAB software where the above optimization algorithms 
are implemented in fmin (1-parameter case) and fmins (n-parameter case) routines. Both algorithms 
involve only calculation of the function to be maximized at certain points but not its derivatives. For the 
multiparameter case the regions of convergence are approximately parallelepipeds. We have summarized 
our results in the Table || below. We have given the values of the intersection of the parallelepipeds with 
the coordinate axes in the parameter space. We have expressed these values in the units of square roots 
of diagonal values of the inverse of the matrix G given by Eq. (|35|) . 
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ro 


r\ 




rz 





10 








1 


0.7 


0.5 






2 


0.2 


0.08 


0.1 




3 


-0.08 


-0.02 


-0.01 


-0.03 



Table 6: Coordinates of the regions of convergence for the polynomial phase signals with s spindowns 
included in the units of the square roots of diagonal elements of the inverse of the matrix G given by Eq. 
(|35|). The region of convergence for the fcth (k = 0, . . . , 3) spindown is the interval [— rjt]. 

In the case of the signal with 3 spindowns our estimation of the radius of convergence is very crude 
because the computational burden to do such a calculation is very heavy. The above results hold for the 
statistics T calculated when data is only signal and no noise. 

In the coarse search we have chosen a rectangular grid in the spindown parameter space with the 
nodes separated by twice the values given in Table ^ and we have chosen the spindown parameter ranges 
to be from —3 to 3 times the square roots of the corresponding diagonal elements of matrix G given by Eq. 
(|35|). We have made 10 4 simulations in the case of a monochromatic signal, 1-spindown, and 2-spindown 
signals and for each signal-to-noise ratio. The case of 3 spindowns turned out to be computationally too 
prohibitive. In each case we have taken the length of the signal to be 2 5 points. 

In our simulations we observe that above a certain signal-to-noise ratio that we shall call the threshold 
signal-to-noise ratio the results of the Monte Carlo simulations agree very well with the calculations of 
the rms errors from the covarince matrix however below the threshold signal-to-noise ratio they differ by 

and has also been observed in 
Eql. There exist more refined 



a large factor. This threshold effect is well-known in signal processing [ 24 
numerical simulations for the case of a coalescing binary chirp signal [ 25 
theoretical bounds on the rms errors that explain this effect and they were also studied in the context of 
the gravitational- wave signal from a coalescing binary p7| . Here we present a simple model that explains 
the deviations from the covariance matrix and reproduces well the results of the Monte Carlo simulations. 
The model makes use of the concept of the elementary cell of the parameter space that we introduced in 
Sec. 3. The calculation given below is a generalization of the calculation of the rms error for the case of 
a monochromatic signal given by Rife and Boorstyn [ p8[ . 

When the values of parameters of the template that correspond to the maximum of the functional 
T fall within the cell in the parameter space where the signal is present the rms error is satisfactorily 
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approximated by the covariance matrix. However sometimes as a result of noise the global maximum is 
in the cell where there is no signal. We then say that an outlier has occurred. In the simplest case we 
can assume that the probability density of the values of the outliers is uniform over the search interval 
of a parameter and then the rms error is given by 

vLt = ^, (132) 

where A is the length of the search interval for a given parameter. The probability that an outlier occurs 
will be the higher the lower the signal-to-noise ratio. Let q be the probability that an outlier occurs. 
Then the total variance a 2 of the estimator of a parameter is the wighted sum of the two errors 

o- 2 = a 2 out q + a 2 CR (l-q) 1 (133) 

where o~cr is the rms errors calculated form the covariance matrix for a given parameter. 

Let us now calculate the probability q. Let T s be the value of T in the cell where the signal is present 
and let T Q be its value in the cells where signal is absent. We have 



/>oo 

1 - q = P{aU: T Q < F s } = / P{all: T Q < T s \T S = T}P{T S = T} dF, 

Jo 



(134) 



where P stands for probability. Since the values of the output of the filter in each cell are independent 
and they have the same probability density function we have 

P{all: T < F S \F S =T}= [P{F < F S \F S = F}}^ 1 , (135) 

where N c is the number of cells of the parameter space. Thus 



1-9=/ Pi{d,T) / p (y)dy 

J Jo 



N c -1 

dT, (136) 



where po and p\ are probability density functions of respectively false alarm and detection given by Eqs. 
@ and (H 



In Figures 10, 11, and 121 we have presented the results of our simulations and we have compared 



them with the rms errors calculated from the covariance matrix. We have also calculated the errors from 



our simple model presented above using Eqs. (|133j ) and ( 136 ). In the case of frequency, spindowns, and 
phase to calculate a out we have assumed uniform probability density. The estimator of the amplitude is 
proportional to the modulus \X\ of the Fourier transform of the data and in the case of the amplitude we 
have calculated a out for the probability density of \X\ assuming that there is no signal in the data. We 
see that the agreement between the simulated and calculated errors is very good. This confirms that our 
simple model is correct. We also give biases of the estimators in our simulations. We see from Figures 
[Io|~[l2] that as signal-to-noise ratio increases the simulated biases tend to zero and the standard deviations 
tend to rms errors calculated from the covariance matrices. 
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A Functions A, B, and C 

The functions A, B, and C in Eqs. (|l4|) for the observation time chosen as an integer number of sidereal 
days take the form (here n is a positive integer) 
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9 cos 4 A cos 4 5 + i sin 2 2A sin 2 26 + ^ (3 - cos 2A) 2 (3 - cos 25) 2 



26 



B 



C 



T a =n2ir/n r 
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-- cos 2 27 (1 + cos 2 A cos 25) 



= 0. 



(137) 

(138) 
(139) 



We see that the functions A, B, and C depend only on one unknown parameter of the signal — the 
declination 5 of the gravitational-wave source. They also depend on the latitude A of the detector's 
location and the orientation 7 of the detector's arms with respect to local geographical directions. 



B The Fisher matrix 



In this appendix we give the explicit analytic formula for the Fisher matrix for the simplified model of 
the gravitational- wave signal from a spinning neutron star. The model is defined by Eqs. ( fill ) an d (^2|) in 
Sec. 4. It has a constant amplitude and its phase is linear in the parameters. In Paper II we have shown 
that this model reproduces well the accuracy of the estimators of the parameters calculated from the full 
model which has amplitude modulation and nonlinear phase. In this paper in Sec. 5 we show that the 
number of templates needed to perform all-sky searches calculated from the linear model reproduces well 
the number of templates calculated from the nonlinear phase model in Ref. |7j- Thus we see that the 
Fisher matrix presented below can be used in the theoretical studies of data analysis of gravitational- wave 
signals from spinning neutron stars instead of a very complex Fisher matrix for the full model. 

In Paper II we have found that the Fisher matrix depends on the choice of the initial time within the 
observational interval (initial time is that instant of time at which the instantaneous frequency and the 
spindown parameters are defined, see Appendix C of Paper II) . However one finds that the determinant 
of the transformation between the two Fisher matrices with different values of the initial time chosen is 
1. Consequently the number of cells and the number of filters do not depend on the choice of initial time. 
We present our analytic formula for the initial time chosen to coincide with the middle of the observation 
interval. This simplifies the analytic expressions considerably. 



The Fisher matrix for all-sky searches with s spindowns included is defined by 
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(140) 



*(i;C) = $ + $(*;£); 



(141) 



the function <I> is given by Eq. (52). We have calculated the Fisher matrix Tf \ for s = 4. The result is 
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The above formulae could further be simplified if we assume that the observation time is an integer 
multiple of one sidereal day. We also note that if we have data corresponding to a full year we can start 
our observation at a time corresponding to any position of the detector in its motion around the Sun. 
This means that in such a case we can choose the phases <f> r and <f> arbitrarily. 

0, . . . , 3 equals to the submatrix of TfX consisting of the first s + 3 



The Fisher matrix rfH for s 



columns and the first s 



3 rows of r^l. 



The reduced matrix rfH defined in Eq. (| 



from the matrix r?'\ by means of the following procedure: take the inverse of T al1 



can also be obtained 
remove the first 



column and the first row of the inverse, take again the inverse of such a submatrix — it equals Tf \ . 

In the case of directed searches the Fisher matrix is also defined by Eq. (14C), but now £ = ($o, £), 



£ = (wq, ...,U s 

included reads 



and the phase $ is given by Eq. ( |78| ) 
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The Fisher matrix T^j with s = 4 spindowns 
\ 
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The Fisher matrix rfe for s = 0, 



, , 3 equals to the submatrix of T d 



(4) 



(143) 



consisting of s + 1 first columns 



and s + 1 first rows of T^j . The reduced matrix rfe defined by Eq. (|6^) can be obtained from the matrix 
rfs) by means of the same procedure as described above for the case of all-sky searches. 



C Suboptimal filtering 

Very often suboptimal filter (or a search template) is proposed in hierarchical signal searches. In such a 
search one passes the data through a suboptimal filter that requires much less computational cost than 
the optimal filter and one registers the condidate events. Then one passes the data through optimal filters 
however only for the values (or around the values) of the parameters of the candidate events to assess 
the significance of the candidate events. In such a search one would like to ensure that there is no loss 
of events. A way to achieve this when using a suboptimal filter is to lower the threshold with respect to 
the threshold chosen for the optimal filter so that the number of expected significant events is the same 
as with the optimum filter. The probability densities derived in Sec. 6.1 can be used to calculate what 
the lowered threshold should be. 

To illustrate the general theory developed in Sec. 6 we have considered the following example. We 
have assumed the observation time T to be 3 days and we have restricted ourselves to directed searches. 
For such a case the model of the phase consistent with the 1/4 of a cycle criterion has s\ = 2 spindowns 
in the dominant term, S2 = 1 spindown in the contribution due to the Earth orbital motion and no 
contribution due to the Earth diurnal motion (s^ — 0), cf. Eq. ( |130| ). We have correlated this signal 
with a template that has s\ = 1, S2 = 1, and S3 = 0. Assuming the gravitational-wave frequency f Q = 1 
kHz and the maximum values of the spindowns for the spindown age r = 40 yr the fitting factor is 0.91, 
the number of cells iV c for the optimal random field is 2.3 x 10 12 and the number of cells N sc for the 
suboptimal random field is 3.7 x 10 12 . We have found that the fitting factor is practically independent 
on the right ascension and the declination of the gravitational-wave source. 

In our computations we assume that we lower the threshold according to the law 

ToL = {To - 1)FF 2 + 1. (144) 

The above rule is motivated by the relation between the expectation value of the statistics T and the 
optimal signal-to- noise ratio given by Eq. (|85|). 

The numerical results obtained using formulae derived in Sec. 6.1 are presented in Figure^. We have 
assumed the false alarm probability to be 1% for the optimal filter. There is one more input parameter 
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that we need in order to calculate the numbers of expected events: the signal-to-noise ratio di for which 
the number of events is 1. In the upper left plot in Figure |l3| we have shown the ratio of the expected 
number of the detected e vents for the suboptimal filtering [calculated from Eq. ( |115| )] and the optimal 
one [calculated from Eq. (114)] as a function of d%. We have assumed that in the suboptimal filter we 
lower the threshold according to Eq. (144). We have also put FF = 0.91. 

To assess the number of events that one loses using a search template Apostolatos p3] assumed that 
the number of detected events decreases as FF 3 . In the right upper plot of Figure 13 we have compared 
the number of detected events calculated from Eq. ( |115| ) and the ones calculated using FF 3 law. We see 
that in general FF 3 law underestimates the event loss. However for the fitting factors close to one the 
difference is small. 

We have calculated the numbers of expected detections and false alarms for the optimal and subopti- 
mal filter both with original and lowered thresholds. The results are presented in the two lower plots in 
Figure [l3| In the plot on the left diamonds mark the ratio of the number of the detected events for the 
suboptimal filter with lowered threshold [calculated from Eqs. (115) and (144)] and the number of events 
detected with the optimum filter [calculated from Eq. (114)]; squares denote the ratio of the number of 
events detected by suboptimal filtering without lowering the threshold and the number of events detected 
with the optimum filter; the solid line gives the fraction of the detected events calculated from FF 3 law; 
all dependencies are shown as functions of the fitting factor. The lower plot on the right gives the ratio of 
the expected number of false alarms with the suboptimal filter and lowered threshold and the expected 
number of false alarms for the optimal filter. 

From our example we see that when using a suboptimal filter by appropriate lowering of the threshold 
we can detect all that events that can be detected with an optimal filter. There is however a limitation to 
threshold lowering arising from the fact that below a certain threshold the false alarm rate can increase 
to an unmanageable level. In the real data analysis there may be other limitations. For example below 
a certain threshold a forest of non-Gaussian events may appear completely obscuring the real signals. 



D The use of parameters a\ and to label the filters 

If one knows the values of the parameters oti, ct2, and f a it is possible to solve Eqs. (|4]) with respect 
to the angles a and S. One can show that each triple [at, ct2, fo) gives two such solutions which can be 
written as follows (note that because S £ [— -|, \ ] to determine S uniquely it is enough to know sin<5): 



sin 5 = p lS ms±Jl-j3f-/3%, (145) 



where 



cos S = VI - sin 2 S, (146) 

Pi — sine sin S ,„ 

sin a = — — , (147) 

cos e cos o 

cos a = (148) 
coso 



Pi ■= 7^, ft:=^. (149) 
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The correspondence between the parameters a\, a2, fo and a, S given by Eqs. (145)-(148) implies 
that one can use a±, 0,2 instead of a, S to label the templates needed for matched filtering. To do this 
the family of templates labelled by a, S (and the other parameters) must be replaced by two template 
families labelled by ax, 012 (and the other parameters). The first family arises when in the original family 



one replaces sin<5, cos5, sina, and cosa by the left-hand sides of Eqs. (145)-(148) with + sign chosen in 



the front of the square root in Eq. (14!:). In the second family the replacements are made with — sign 
chosen. The filters labelled by parameters a,\ and «2 will to a good appproximation be linear and the 
theory of data processing developed in this paper applies to such a filtering scheme. 

When as a result of filtering of the data one gets a significant event one obtains at the same time 
the maximum likelihood estimators of the parameters a\, «2, fo (and the others). One can obtain the 
maximum likelihood estimators of the position (a, S) of the gravitational- wave source in the sky by means 
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of Eqs. ( 145 )— ( |14q ). Note that one should expect to get the maximum correlation for a template belonging 
to one out of two families described above, what means that after filtering one would also know which 
sign on the left-hand side of Eq. fll45| ) should be chosen. 

The covariance matrix for the parameters a, 6, and f a can be obtained from the covariance matrix 
for the parameters a%, a2, and f Q by means of the law of propagation of errors. Let us introduce 



x := (ai,a 2 ,f ), y := (a,<$, / ). 



(150) 



Let C x be the covariance matrix for the parameters x, then the covariance matrix C y for the parameters 
y can be calculated as follows: 



(151) 



where T denotes matrix transposition and the Jacobi matrix J has components: 
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All derivatives entering Eq. ( |152| ) can be calulated using Eqs. (145)— (148). 
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Figure 1: Probability density of the false alarm (upper plot) and the false alarm probability (middle 
plot) of 127 bins of the statistics Tu (for k = 2,..., 128) given by Eq. (p0[). The false alarm probability 
of the first (k — 1), zero frequency bin and the last (k = 129), Nyquist frequency bin is given in the left 
lower and the right lower plot, respectively. The continuous lines in the upper and the middle plots are 
theoretical distributions given by Eqs. ( p7j ) and (^9|) for n = 2 and the continuous lines in the lower plots 
are theoretical distributions that follow from the cumulative % 2 distribution with 1 degree of freedom. 
In the simulation one million of sequencies of 256 random independent samples drawn from zero mean 
and unit variance normal distribution were generated and modulus of their discrete Fourier transforms 
evaluated. The results of the simulations are marked by the squares. 
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Figure 2: False alarm probabilities for a monochromatic (three upper plots) and a linearly frequency 
modulated signals (three lower plots). The same random sequences of length N = 2 s were generated 
as in the simulation in Figure 1 except that experiment was repeated 10 5 times. The results of the 
simulation are marked by the squares (no zero padding) and by the circles (for the signal padded with 
37V zeros). The ratio Tp (cells) is the quotient of the false alarm probability obtained from the simulations 
and calculated from Eq. (^8|) whereas rp(Euler) is the quotient of the false alarm probability obtained 
from the simulations and calculated from Eq. (Ejj). 
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Figure 3: Probability of detection (plots on the left) and the receiver operating characteristic (plots on 
the right) for a monochromatic (upper plots) and a linearly frequency modulated signal (lower plots). 
The same random sequences of length N = 2 8 were generated as in the simulation in Figure 1 except 
that the experiment was repeated 10 times. The results of the simulation are marked by the circles. 
Theoretical distributions are given by solid lines. Probability of detection is calculated from Eqs. ( p8| ) 
and ( |30| ) for n = 2 and optimal signal-to-noise ratio d — 4. The receiver operating characteristics are 
parametric curves with signal-to-noise ratio d as a parameter, they are calculated from Eqs. (|3(]) and ( ^8|) 
for d = 4, 5, and 6. 
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Figure 4: Volume of one cell (in units s -2 ) in all-sky searches as a function of the observation time T 
for the LIGO Hanford detector (latitude A = 46.45°). We have calculated the volume of one cell from 
Eq. ([75|). The lines shown in the plot correspond to different numbers s of spindowns included: s = 4 
(solid), s — 3 (dotted), s — 2 (dashed), s — 1 (dotted/dashed), and s = (double dotted/dashed). We 
have set cj) r = 1.456 and <f) = 0.123. 
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Figure 5: Number of cells in all-sky searches as a function of the observation time T Q for different values 
of the minimum spindown age T min and the maximum gravitational-wave frequency / max (the minimum 
gravitational- wave frequency / m i n = 0). The lines shown in the plots correspond to different numbers s 
of spindowns included: s = 4 (solid), s = 3 (dotted), s = 2 (dashed), s = 1 (dotted/dashed), and s = 
(double dotted/dashed). We have assumed the LIGO Hanford detector and we have put (p r — 1.456 and 
(f> = 0.123. 
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Figure 6: Dependence of number of cells in all-sky searches on the angles <f> r , <j> , and the latitude 
A of the detector's site. We have chosen the minimum spindown age T m i n = 40 years, the maximum 
gravitational-wave frequency / max = 1 kHz, and the minimum gravitational-wave frequency / m ; n = 0. 
The plots (a), (c), and (e) are for the observation time T — 7 days (and the number of spindowns s = 2); 
(b), (d), and (f) are for T = 120 days (and the number of spindowns s = 3). In the plots (a), (b), (c), 
(d) we have used the latitude A = 46.45° of the LIGO Hanford detector; in (a), (b), (e), (f) we have put 
(f> — 0.123; and in (c), (d), (e), (f) we have used 4> r — 1.456. 
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Figure 7: Number of cells in directed searches as a function of the observation time T for different values 
of the minimum spindown age r m i n and the maximum gravitational-wave frequency /max (the minimum 
gravitational- wave frequency / m in = 0). The lines shown in the plots correspond to different numbers s 
of spindowns included: s = 4 (solid), s = 3 (dotted), s = 2 (dashed), s = 1 (dotted/dashed), and s = 
(double dotted/dashed). 
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Figure 8: Number of filters in all-sky searches as a function of the observation time T Q for different 
values of the minimum spindown age r m i n and the maximum gravitational- wave frequency / max - The 
lines shown in the plots correspond to different numbers s of spindowns included: s = 4 (solid), s = 3 
(dotted), s = 2 (dashed), s = 1 (dotted/dashed), and s = (double dotted/dashed). We have assumed 
the LIGO Hanford detector and we have put <j) r = 1.456 and <p = 0.123. 
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Figure 9: Number of niters in directed searches as a function of the observation time T for different 
values of the minimum spindown age r m i n and the maximum gravitational- wave frequency /max- The 
lines shown in the plots correspond to different numbers s of spindowns included: s = 4 (solid), s = 3 
(dotted), s = 2 (dashed), s = 1 (dotted/dashed). 
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Figure 10: Simulations of the biases (plots on the left) and the rms errors (plots on the right) for a 
monochromatic signal. The results of the simulations are marked by the circles. The £-axes are labelled 
by the optimal signal-to-noise ratio. The thin solid line s in t he plo ts o n the right are calculated from the 
covariance matrix and the thick lines follow from Eqs. (133) and (|l3(j). 
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Figure 11: Simulations of the biases (plots on the left) and the rms errors (plots on the right) for a 
1-spindown signal. The results of the simulations are marked by the circles. The x-axes are labelled by 
the optimal signal-to-noise ratio. The thin solid lines in the plots on the right are calculated from the 
covariance matrix and the thick lines follow from Eqs. (133) and (p~3q). 
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Figure 12: Simulations of the biases (plots on the left) and the rms errors (plots on the right) for a 
2-spindown signal. The results of the simulations are marked by the circles. The x-axes are labelled by 
the optimal signal-to-noise ratio. The thin solid lines i n th e plot s on the right are calculated from the 
covariance matrix and the thick lines follow from Eqs. (|l33| ) and (|l3(j ). 
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Figure 13: Suboptimal filtering. In the upper left plot we show the ratio N ^((iisub, Fol) / 'Noidi, T ) of 
the expected number of the detected events for the suboptimal filtering [calculated from Eq. (115)] and 
the optimal one [calculated from Eq. (114)] as a function of the signal-to- noise ratio d\ (the signal-to- noise 
ratio for which the number of events is one). We have assumed that in the suboptimal filter we lower the 
threshold according to Eq. (144). We have also put FF = 0.91. In the right upper plot we give the ratio 
t = [N S £)(di su b, F l) /Nd(cIi, Fo^/FF 3 as a function of the fitting factor (we have used d\ — 16.6). In 
the left lower plot diamonds mark the ratio N s i>(di su b, J- l) /N[>(di, J- ) of the number of the detected 
events for the suboptimal filter with lowered threshold [calculated from Eqs. (115) and (144)] and the 
number of events detected with the optimum filter [calculated from Eq. (114)]; squares denote the ratio 
N s d (rfisub; 3 r )/N]j(di, T ) of the number of events detected by suboptimal filtering without lowering the 
threshold and the number of events detected with the optimum filter; the solid line gives the fraction of 
the detected events calculated from FF 3 law; all dependencies are shown as functions of the fitting factor 
(we have put d\ = 16.6). The lower plot on the right gives the ratio of the expected number of false 
alarms with the suboptimal filter and lowered threshold and the expected number of false alarms for the 
optimal filter. 
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